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FOREWORD 


The present report is one of a series of six reports, published simul- 
taneously, which describe analyses and computational procedures for: 1) pre- 

diction of the in-depth response of charring ablation materials, based on one- 
dimensional thermal streamtubes of arbitrary cross-section and considering 
general surface chemical and energy balances, and 2) nonsimilar solution of 
chemically reacting laminar boundary layers, with an approximate formulation 
for unequal diffusion and thermal diffusion coefficients for all species and 
with a general approach to the thermochemical solution of mixed equilibrium- 
nonequilibrium, homogeneous or heterogeneous systems. Part I serves as a 
summary report and describes a procedure for coupling the charring ablator 
and boundary layer routines. The charring ablator procedure is described in 
Part II, whereas the fluid-mechanical aspects of the boundary layer and the 
boundary- layer solution procedure are treated in Part III. The approximations 
for multicomponent transport properties and the chemical state models are 
described in Parts IV and V, respectively. Finally, in Part VI an analysis • 
is presented for the in-depth response of charring materials taking into ac- 
count char-density buildup near the surface due to coking reactions in depth. 

The titles in the series are: 


Part I 

Part II 

Part III 
Part IV 

Part V 

Part VI 


Summary Report: An Analysis of the Coupled Chemically Reacting 

Boundary Layer and Charring Ablator, by R. M. Kendall, E. P. 
Bartlett, R. A. Rindal, and C. B. Moyer. 

Finite Difference Solution for the In-depth Response of Charring 
Materials Considering Surface Chemical and Energy Balances, by 
C. B. Moyer and R. A. Rindal. 

Nonsimilar Solution of the Multicomponent Laminar Boundary Layer 
by an Integral Matrix Method, by E. P. Bartlett and R. M. Kendall. 

A Unified Approximation for Mixture Transport Properties for Multi- 
component Boundary-Layer Applications, by E. P. Bartlett, R. M. 
Kendall, and R. A. Rindal. 

A General Approach to the Thermochemical Solution of Mixed Equilib- 
rium-Nonequilibrium, Homogeneous or Heterogeneous Systems, by 
R. M. Kendall. 

An Approach for Characterizing Charring Ablator Response with In- 
depth Coking Reactions, by R. A. Rindal. 


This effort was conducted for the Structures and Mechanics Division of 
the Manned Spacecraft Center, National Aeronautics and Space Administration 
under Contract No. NAS9-4599 to Vidya Division of Itek Corporation with Mr. 
Donald M. Curry and Mr. George Strouhal as the NASA Technical Monitors. The 
work was initiated by the present authors while at Vidya and was completed 
by Aerotherm Corporation under subcontract to Vidya (P.0. 8471 V9002) after 
Aerotherm purchased the physical assets of the Vidya Thermodynamics Depart- 
ment. Dr. Robert M. Kendall of Aerotherm was the Program Manager and Prin- 
cipal Investigator. 
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ABSTRACT 


A laminar nonsimilar boundary- layer procedure is described which yields 
accurate solutions for a broad range of problems, in its current formula- 
tion, solutions can be obtained for any equilibrium chemical environment with 
specified rate-controlled reactions at the surface. It has been used to treat 
a variety of ablating and nonablating surface boundary conditions including 
coupled energy and mass balances. The formulation considers unequal diffu- 
sion and thermal diffusion coefficients for all species in a particularly 
convenient manner through a bifurcation approximation for binary diffusion 
coefficients. The multicomponent viscosity and thermal conductivity of the 
mixture are determined by use of Sutherland-Wassiljewa type approximations. 

The procedure is readily applicable to inclusion of one-dimensional radia- 
tion emission and absorption and a general nonequilibrium chemical model. 

The procedure combines features of the general integral relations approach 
with those of matrix solution techniques. Following the former, smooth func- 
tions (in particular, cubic spline functions) are chosen to relate the princi- 
pal dependent variables to their derivatives. This enables the attainment of 
an accurate solution with relatively few entries into the conservation equa- 
tions (3 to 4 place accuracy with 7 to 11 spline points) . From the latter, 
the concept of treating the entire solution as a set of simultaneous, non- 
linear algebraic equations is adopted. This technique results in linearized 
coupling between all relations required to characterize the boundary layer, 
and thus assures a general, rapid, and stable iterative convergence. As a 
consequence, the damping of corrections has seldom been required. Computa- 
tional speed appears to be an attractive feature based on the few comparisons 
with other techniques which have been possible to date. 
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AN INTEGRAL-MATRIX METHOD FOR NONSIMILAR 
SOLUTION OF THE MULTICOMPONENT LAMINAR BOUNDARY LAYER 

SECTION 1 
INTRODUCTION 

A computational procedure is described which is suitable for obtaining 
accurate numerical solutions of the nonsimilar multicomponent laminar bound- 
ary layer with arbitrary equilibrium or nonequilibrium chemical systems, 
unequal diffusion and thermal diffusion coefficients for all species, radia- 
tion absorption and emission, and a variety of surface boundary conditions 
including intimate coupling with transient charring-ablation energy and mass 
balances . A Fortran IV computer program has been developed in accordance 
with this analysis with the exceptions that 1) the chemical system is presently 
limited to equilibrium, with or without selected rate-controlled surface re- 
actions or surface catalyzed reactions, and 2) radiation absorption and emis- 
sion is not currently permitted. This computer program, designated BLIMP, 
for Boundary Layer Integral Matrix Procedure, is described in Ref. 1. 

The computational procedure has been developed while attempting to take 
advantage of the most attractive features of other boundary-layer procedures . 

In light of the application of the procedure to be adopted, certain specific 
requirements seemed appropriate. In particular, minimization of the number 
of “nodal points" required to obtain a solution was judged to be of prime 
importance as a consequence of the relatively large times associated with 
state calculations for a general chemical environment and, in the streamwise 
direction, because of the desire to couple the boundary layer procedure to a 
transient internal conduction or ablation solution. 

For a given accuracy, the number of necessary “nodal points" in the sur- 
face normal direction is controlled primarily by the nature of the functions 
which relate the dependent variables (and their derivatives) to the indepen- 
dent variable. Thus the continuous functions typically used in integral 
relations approaches require fewer "nodal points"* than the discontinuous 
functions implied by most finite difference approximations, in order to per- 
mit relatively flexible profiles, sets of connected cubics were selected to 
represent enthalpy, velocity, and elemental concentrations. The first and 
second derivatives of these cubics were made continuous at the connecting 
points. The advantages of such a "spline fit" are considered, for example, 
in Ref. 2. 


*The term "nodal point" is meant to encompass the integral strips of Pallone 
and the matching points used by Dorodnitsyn. 




If the general integral relations approach is followed, weighting func- 
tions must be selected. In the present study this selection was based pri- 
marily on the complexity of the resultant algebra. Studies were made using 
Dirac delta weighting functions (i.e., a differential approach*) and step 
weighting functions similar to those used by Pallone 3 which indicated,** when 
other aspects of the procedure were unchanged, no definite superiority in 
terms of accuracy or stability. Because all of the complexities introduced 
by the generalization of the thermodynamic and transport properties of the 
system occur within a divergence term, step weighting functions produce 
markedly simpler algebra and, consequently, were adopted for the present 
procedure . 

In the past when relatively large spacing in the streamwise direction 

has been desired, iterative procedures have generally been used to assure 
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accuracy and stability. In many instances ' these procedures have treated 
the solution in a manner resembling those used for similar solutions but with 
the addition of finite difference representations for the nonsimilar terms, 
a procedure which eliminates the necessity of special starting techniques. 
Using this basic approach, the specific treatment adopted in the current 
study follows most closely the matrix procedure used by Leigh 6 wherein the 
iteration is a consequence of the solution of a set of linear and nonlinear 
algebraic relations. Whereas a special successive approximation procedure 
was used by Leigh, the general Newton-Raphson technique was used in the pres- 
ent procedure. This technique results in linearized coupling between all 
relations required to characterize the boundary layer, and thus assures a 
more general, rapid and stable iterative convergence. 

The present document concentrates on the fluid mechanical aspects of 
the problem and describes the basic numerical solution procedure. The pro- 
cedures employed for calculating the equilibrium state of the gas and sug- 
gested for including rate-controlled reactions are described elsewhere 7 
since they are conveniently treated as subroutines to the basic boundary 
layer computational procedure. However, the terms which are directly in- 
volved in the boundary layer equations such as the "elemental source term" 
which arises from kinetic considerations are included in the present develop- 
ment. Similarly, radiation absorption and emission enters directly into the 
conservation equations only as a net radiation flux term in the energy equa- 
tion. The calculation of this term could also be conveniently accomplished 
by a subroutine. A one-dimensional model for net radiation flux which rep- 

Q 

resents an extension of the work of Cess to allow an angular-dependent in- 
cident radiation flux at the boundary- layer edge is presented in Appendix E. 

4 

*This correspondence is pointed out by Dorodnitsyn. 

**The results of these studies are discussed in an appendix to this report. 
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Multicomponent transport properties are based on a newly developed approxima- 
tion described in Ref. 9. Modification of the conservation equations as a con- 
sequence of this approximation is described herein. Finally, the procedures 
employed for coupling to a transient charring ablation program are described 
in Ref. 10. 

The governing differential equations for laminar or turbulent flow are 
presented in Section 2. The laminar form of the equations are normalized by 
a modified Levy-Lees transformation in Section 3. The modification consists 
of a coordinate stretching parameter which permits the establishment of an 
efficient universal boundary layer nodal network. In Section 4, the trans- 
formed conservation equations are integrated and the connected-cubic func- 
tional relationships are introduced through truncated Taylor series expan- 
sions. The procedure utilized to solve these equations is described in 
Section 5. First, the Newton-Raphson linear recurrence formulas are devel- 
oped. A matrix reduction procedure is then described which takes full 
advantage of the linear Taylor series expansions and simplifies the gener- 
alization of surface boundary conditions. 

In Section 6, comparisons to other numerical solutions are shown for 
several uncoupled nonreacting boundary layer problems. Generally, 3- to 4- 
place accuracy is obtained for 7-point boundary-layer solutions, and the 
solution usually converges in 3 or 4 iterations. Solutions for chemically 
reacting boundary layers are presented in Section 7. 

SECTION 2 

BOUNDARY LAYER CONSERVATION EQUATIONS 

In this section are presented the differential equations which govern 
laminar or turbulent flow in a planar or axisymmetric compressible boundary 
layer with mass addition, equilibrium or nonequilibrium chemical reactions, 
multicomponent diffusion, thermal diffusion, and radiation. Unequal diffu- 
sion and thermal diffusion coefficients for all diffusing pairs are in 
accordance with the unified approximation presented in Ref. 9. The equations 
derived are essentially an extension of those derived in Ref. 11 to include 
radiation, thermal diffusion, and unequal binary-diffusion coefficients. The 
diffusion introduced by pressure gradients and body forces are neglected. 

The standard definitions of time-averaged turbulent quantities and rela- 
tive order of magnitude are employed (Refs. 11 and 12) . The turbulent trans- 
port terms are expressed in the Boussinesq form, that is, eddy viscosity, 
eddy diffusion, and eddy conductivity. Hence, all the terms in the equations 
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are time-averaged quantities and no need exists for using a superscript bar.* 
In the order-of-magnitude arguments, terms of the following types have been 
eliminated: 1) triple correlations, 2) derivatives of turbulent correla- 

tions parallel to the wall, and 3) correlations involving turbulent compo- 
nents of molecular transport mechanisms. 

A mass balance of an individual species in a unit volume results in the 
relationship 


K ds 
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where s and y are the streamwise and normal coordinates, respectively, 
u and v are the velocity components in the s and y directions, respec- 
tively, is the mass fraction of species i, r Q is the radius of the 

body in a meridian plane for an axisymmetric shape, k is zero for a flat 
plate and unity for a body of revolution, p is the density, represents 

the rate of mass generation of species i per unit volume due to chemical 
reaction, pe Di is defined in terms of the correlation of the fluctuating 
components of concentration and normal velocity, that is. 


(pv)' k 7 

pe D ± = SKT7§y 
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and 


j_^ is the mass-diffusion rate of species i due to molecular processes. 
When Eq. (1) is summed over all the species in the system, utilizing 


i i 

which results from the definition of mass diffusion, and utilizing 
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*Accordingly, pv represents pv, not p v. 
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which results from conservation of mass, there results 
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which is the familiar global continuity equation. 

When Eq. (3) is considered together with Eq. (1) , the more conventional 
species conservation equation is obtained: 
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The streamwise momentum equation can be written as 
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where P is the pressure and the eddy viscosity is defined in terms of the 
Reynolds stresses of turbulent flow by 


pe 
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The momentum equation for forces and fluxes normal to the surface is 
given by 


3P\ _ £ul 

dy 1 = - 


(7) 


where r is the radius of curvature of a surface streamline, 
c 

The energy equation for this general system is 
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(equation continued on next page) 


5 



(equation continued from previous page) 



where H,p is the total enthalpy (static plus kinetic) 
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h is the static enthalpy including chemical as well as sensible contribu- 
tions 
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h^ is the static enthalpy of species i 
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T is the temperature, h? 
the specific heat of species 
gaseous mixture 


is the heat of formation of species 

i* 
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C is the frozen specific heat of the 
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X is the thermal conductivity, R is the gas constant, is the mole 

fraction of species j, 7%. is the molecular weight of species i, J}. . is 
the binary diffusion coefficient of species i into j, is the multi- 

component thermal diffusion coefficient of species i, the turbulent enthalpy- 
transport coefficient is defined by 


Y vp^’v 
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and q r is the net one-dimensional energy flux towards the surface due to 
radiation absorption and emission.* In Eq. (8) the turbulent contribution 
to the DuFour effect (the double summation term) has been neglected since 
significant diffusion thermo occurs only in the laminar region where tempera- 
ture gradients are severe. 

When the assumption of equal diffusion coefficients is made, a substantial 
simplification of the problem results if the species conservation Equations 
(4) are multiplied by a^, defined as the mass fraction of element k in 
species i, and the resulting terms are summed over all species (known 
as the Shvab-Zeldovich transformation) . When this is done, there is a re- 
duction in the number of conservational equations from the number of species 
(typically 20 to 50) to the number of elements (usually 2 to 6) . In addition, 
the resulting equations are simplified since the source terms are eliminated. 
Furthermore, elements vary more smoothly across the boundary layer than do 
the molecular species, and hence are better represented numerically. To illus- 
trate, when all binary-diffusion coefficients are assumed equal and in the 
absence of thermal diffusion, the j i can be expressed by Pick's law 

SK. 

j i = " pi> i2 aF (14) 


Substituting this into Eq. (4) and performing the Shvab-Zeldovich transforma- 
tion results in the following elemental conservation equations for the laminar 
or turbulent boundary layer: 
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where is the mass fraction of element k in the system defined by 

\ s Y (16) 

i 

It has also been assumed that all e D = e D . 

When diffusion coefficients are not equal, Fick's law does not apply. 

The diffusional fluxes, j., must then be expressed in terms of multicomponent 
diffusion coefficients, 

*A model for q r which allows an angular-dependent incident radiation flux at 
the boundary-layer edge is developed in Appendix E. 
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or via the Stefan-Maxwell relations 
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Utilization of the Stef an- Maxwell equations in conjunction with the species 
conservation equations is awkward even in the absence of thermal diffusion 
effects, since the diffusional flux, j^, is expressed implicitly in terms of 
mole fractions and their gradients. Hence, use is often made of Eq. (17) 
together with the multicomponent diffusion coefficients, for example, in 
Refs. 14 through 16. However, each of the (I s -I) multicomponent diffusion 
coefficients depends upon local concentrations and upon (I 2 -l)/2 symmetric 
binary diffusion coefficients, # where X is the total number of species 
being considered. 

A bifurcation approximation to binary diffusion coefficients introduced 
17 

by Bird and utilized herein permits explicit solution of the Stef an-Maxwell 
relations for j ^ in terms of gradients and properties of species i and 
of the system as a whole. The approximation can be expressed in the form: 

« D/P ± F j (19) 

with D(T,p) a property of the given multicomponent mixture and F. (T) a 
property of the species in the mixture.* It is apparent when consider- 

ing more than 3 species** that Eq. (19) is indeed approximate, since (I 2 -l)/2 
diffusion coefficients, .5—# are replaced by I diffusion factors, F^. 
Equation (19) should thus be viewed as a correlation equation for actual 
binary diffusion coefficient data. The are determined for a given chem- 

ical system by a least-squares fit of actual diffusion data. 


*The effect of pressure can be absorbed entirely into the D since .&• . is 
inversely proportional to pressure. It will be shown that the dominant; 
temperature effect can also be absorbed into the D so that the F^ are 
nearly constants for a given molecular set. 

**Eq. (19) is exact for a ternary system. 
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The accuracy of the correlation was investigated by Bird^ 7 for a five 
component mixture containing hydrogen and shown to be surprisingly good, the 
maximum error in any being 4 percent. In order to establish more gen- 

erally the adequacy of the approximation, correlations were performed for 

g 

several chemical systems including a 16-component (120 j^j) C-H-O-N system. 
These studies have demonstrated that Eq. (10) has general applicability, and 
represents the ih ^ for nearly all diffusing pairs within 5 percent. The 
largest single error in obtained in these correlations has never ex- 

ceeded 15 percent or so. 

Introducing this approximation into the Stef an- Maxwell relations, it is 
shown in Ref. 9 that the j . can be expressed explicitly as 
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where is a quantity which for unequal diffusion lies between a mass and 

a mole fraction and is defined by 
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v 2 


( 21 ) 


and and ^ are system quantities defined by 



( 22 ) 


It can be seen from Eqs . (21) and (22) that £ Z ± = 1. When diffusion coef- 
ficients are assumed equal, setting = 1 yields = K^, m = 1, and H2 - 

In Ref. 9 it was observed that the are weak functions of tempera- 

ture. Thus Eq. (20) can often be simplified to 
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Substituting Eq. (23) into the boundary-layer species conservation equa- 

q 

tion (Eq. (4)) and performing the Shvab-Zeldovich transformation yields the 
following conservation equations for chemical elements in a multicomponent, 
laminar or turbulent boundary layer with unequal diffusion coefficients: 
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where is defined by Eq. (16) and 
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The term 0^ requires some discussion. Introduction of the Shvab-Zeldovich 
transformation eliminates the chemical production terms i|f when the boundary 
layer is everywhere in local equilibrium since 0^ is then equal to zero. With 
the introduction of nonequilibrium, the approach of Ref. 7 generalizes the term 
"element" and results in an expanded (in terms of "elements", k) array, tak- 

ing advantage of all equilibrium aspects of the system. As a consequence, the 
0^ may no longer equal zero and thus cannot, in general, be omitted from Eq. 

(24) . This general mixed equilibrium, nonequilibrium approach results in more 
equations of the form of Eq. (24) than in the purely equilibrium system. Except 
in the limit of all reactions being kinetically controlled, the number of equa- 
tions of the form of Eq. (24) is, however, less than the number of molecular 
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conservation equations (Eq. (4)). The local state of the gas for this nonequi 
librium system is defined in terms of the expanded set of "elements 0 . The pro 
duction or destruction rates, 0^, of these "elements" becomes a state property 
and can be evaluated along with other local system properties. 

It should be noted that the Shvab-Zeldovich transformation is still 
possible without the approximation for embodied in Eq. (19) , but solu- 

tion then depends upon (I s -I) multicomponent diffusion coefficients, each 
of which depends upon (I s -I) /2 symmetric binary diffusion coefficients and 
upon concentrations of all species. Therefore, use of the approximation for 
iK j embodied in Eq. (19) should be looked upon as a computational convenience. 

At this point the multicomponent thermal diffusion coefficients, D^ T , 

still appear in the conservation relations. Theoretical equations for D. T 

13 1 

are quite complicated and these have to be solved at every boundary-layer 

point since they are strongly concentration dependent. Therefore, a correla- 
tion of binary thermal diffusion data was conducted which yielded, upon gen- 
eralization to multicomponent systems, the following simple relation: 
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T 

with the empirical constant c, about -0.5. This approximation for D. 

. . r 13 1 

satisfies the requirement that they sum to zero, the observation that they 

are independent of fluxes, and the assumption that thermal diffusion of spe- 
cies i should behave nearly as though it were in a system of species i 
and a species representative of the mixture as a whole. The approximation 
represents binary thermal diffusion data reasonably well (within 10 percent 
or so, considering a wide range of molecular weights and variation of mass 

fractions from zero to 100 percent) . An accuracy study of multicomponent 
T 

has not yet been accomplished. However, the generalizations which were 

employed appear to be in basic accord with the approximate model for multi- 

18 

component thermal diffusion coefficients developed by Laranjeira. 

Inserting this approximation into Eqs. (24) and (8) and performing all 
summations yields the following relations for diffusive mass flux of species 
i, j i# diffusive mass flux of element k, j^, and diffusive heat flux, q a , 
respectively 
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The elemental species conservation equation thus becomes 
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while the energy equation can be expressed as 
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Equations (27a) through (27c) are derived in Appendix A. For assumed equal 
diffusion coefficients, = 1/7/1, = C and K = h. When thermal diffu- 
sion is to be neglected - 0 and p 4 = Jtn |J 2 . 

Equations (3), (5), (7), (29), and (30) comprise the boundary-layer 

conservation equations incorporating the approximations for unequal thermal 
and multicomponent diffusion coefficients embodied in Eqs . (19) and (26). It 
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should be emphasized that the numerical solution procedure described in a 
later section is not dependent upon the use of these convenient approxima- 
tions . 

Any consistent set of boundary conditions can be applied which yield, in 
effect, u/u^, an< 3 the and their first derivatives at the edge of the 

boundary layer, and wall values of u/u , mass flux, H_ (or its gradient), and 
the (or their gradients) . Specific boundary conditions will be introduced 

in Section 3 . 

The mathematical specification of the boundary- layer problem is completed 
by introduction of the remaining multicomponent transport properties, the equa- 
tion of state for a gaseous mixture, the equilibrium relations, and the kinetic 
relations. These are described in Refs. 7 and 9. 


SECTION 3 

THE TRANSFORMED NON$IMILAR LAMINAR BOUNDARY- LAYER EQUATIONS 

From the original formulations of Blasius, a continuing effort has been 
expended in the search for more general means of reducing the partial-differ- 
ential equations of the boundary layer to total-differential equations. Bas- 
ically, this involves the search for a new coordinate system (r|,§) related 
to the original system (y,s) and certain of the dependent variables, in which 
the §-wise variations of functions of the dependent variable either vanish 
or become of second order. A successful similarity transformation, as this 
is called, results in ^-derivatives vanishing, but this occurs only under 
certain conditions which are generally quite restrictive. Currently, the 
most popular transformation represents a combination of the Levy and Mangier 
and the Howarth-Dorodnitsyn transformations. This particular form was sug- 
gested by Lees (Ref. 19) among others, and is known by a variety of names 
including Lees-Dorodnitsyn, Levy-Lees, Mangler-Dorodnitsyn, and Dorodnitsyn- 
Stepanov. This transformation is as follows: 
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where, in this and subsequent equations, the subscripts w and e refer, 
respectively, to the wall and to a reference condition which can be taken as 
the boundary-layer edge in the absence of an entropy layer (to be discussed) . 
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In this section, a slightly modified form of this transformation is 
applied to the laminar form of the boundary-layer conservational equations 
presented in Section 2. Although the boundary- layer equations remain partial- 
differential equations when the nonsimilar terms are retained, the Levy-Lees 
transformation is still quite advantageous, since it aids in the specifica- 
tion of the boundary conditions, it eliminates the global conservation equa- 
tion from the set of relations to be solved, and it normalizes the boundary- 
layer thickness. In addition, the nonsimilar terms are often small; hence, 
they can be investigated individually and eliminated for certain classes of 
problems . 

If the conventional Levy-Lees transformation embodied in Eqs . (31) and 

(32) is utilized, the transformed boundary- layer thickness is uniform for a 
similar boundary layer. However, when the boundary layer is highly non- 
similar (e.g., as a result of large blowing or suction, severe pressure gradi- 
ents, or surface discontinuities) the transformed boundary-layer thickness 
can vary by a factor of two or more. 'Therefore, it is useful to normalize 
the boundary layer further by stretching the r\ coordinate: 


§ = ? 



(33) 


where is a function of § only and is determined implicitly during the 

numerical solution. This makes possible the efficient use of a universally 
applicable nodal network which can be chosen a priori once and for all. The 
use of such a universal nodal network is highly desirable as the linearity of 
a large body of equations (Taylor series expansions of primary variables) 
during the numerical solution procedure is retained. In addition, it reduces 
the variation of boundary-layer parameters along a grid line from one stream- 
wise station to the next. 

Since a new variable <* H (T) is introduced, an additional relation is re- 
quired. This is conveniently supplied by constraining an internal nodal point 
near the boundary-layer edge, T} c , to have a specified streamwise velocity, 
c, near (but something less than) the edge value: 
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(34) 
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where f is the stream function defined as 


f - 




( 35 ) 


and the prime denotes partial differentiation with respect to rj, so that 


f 1 



(36) 


To illustrate, selection of r\ = 4.0, r\ c = 2.4, and c = 0.90 yields a R « 1.0 
for the Blasius problem (incompressible flow along a flat plate at zero in- 
cidence with no mass addition). The f 1 | I s utilized in Eq. (34) in anti- 
cipation that the u / u e may not be unity at the edge of the boundary layer 
for superorbital reentry problems involving an entropy layer or nonadiabatic 
flow field. 

In order to illustrate further the procedure, two extreme velocity dis- 
tributions are compared in Figure 1 to the Blasius solution. The profiles 
shown are those for a boundary layer near separation and one near blowoff. 

These were calculated using the integral matrix method for numerical solution 
of the boundary layer (described later) . It can be seen that by the use of 
the transformation (Eq. (33)) together with the arbitrarily chosen constraint 
(Eq. (34)) that the boundary-layer edge occurs at about the same value of r\ 
for the three problems. It should be noted that this is accomplished with 
little mathematical complexity? only two terms involving derivatives of a R 

appear (both in the momentum equation) . 

Variation of a„ with Reynolds number (proportional to §) is shown in 
ri 

Figure 2 for a nonsimilar boundary layer with constant blowing and one with 
constant suction. These results were also obtained using the integral matrix 
technique. The behavior of a u is indicative of the increase (and decrease) 
of the T| at the edge of the boundary layer. The desirability of the stretch- 
ing transformation is made apparent by this example. Without the use of the 
transformation to r\, the choice of an ri v sufficiently large to charac- 
terize accurately the boundary layer at large distances from the leading edge 
would be inefficient near the leading edge. Furthermore, it would be re- 
quired to make an estimate for r l max as It is not known a priori. 
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Figure 2. Variation of ajj with Reynolds number for Incompressible Boundary 
Layers with Constant Blowing and Suction, p v /p u = 1 x 10 -3 






Transformation of the independent variables s and y into the Levy- 
Lees variables § and r\ is conveniently accomplished through the use of 
the operator 
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In addition, the partial derivatives are given by 
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Equations (37) through (41) are derived in Appendix B. 

Utilization of Eqs. (37) through (41) results in the following trans- 
formed equations for diffusive fluxes and conservation of momentum, elemental 
species, and energy in a laminar compressible nonsimilar boundary layer with 
mass addition, chemical reactions, and unequal diffusion and thermal diffu- 
sion coefficients (using the approximations embodied in Eqs. (19) and (26)). 
These equations are derived in Appendix B. Throughout the remainder of 
this document the bar is dropped from 5 # h and the prime refers to partial 
differentiation with respect to r\ except when noted otherwise: 

■hVi 

Diffusive flux of i species 


j± = a.*j i * 

where j^* is a normalized diffusive flux of species i 
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a* is the normalizing parameter defined by 
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C is defined by 
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and Sc is a system property defined by 
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The Sc is a Schmidt number based on the self diffusion coefficient for a 
ficticious species representative of the system as a whole. 


Diffusive flux of elemental species 


j k = a*3y* 


where j k * is a normalized diffusive flux of element k 
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Diffusive heat flux 
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where q * is a normalized diffusive heat flux 
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1 C P + u 1 m 2 


+ c t RTH 3 ' + (h - h + c t RTU 3 )u 4 ' 


(50) 
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and Pr is the Prandtl number based on the frozen specific heat 


Pr 



X 


Streamwise momentum equation 


ff " + 
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• 2?a H 3 1 
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|2 H 
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where p is the streamwise pressure-gradient parameter 

d in u_ 


3 = 2 


d in ? 


( 51 ) 


(52) 


(53) 


Normal momentum equation 


P* 


u e ( 2?) 1 / 3 (f’) s 


a H r c r o 


0 


(54) 


In the present study it will be assumed that the pressure is constant across 
the boundary layer. Equation (54) is therefore replaced by that statement, 
and the partial derivatives of pressure in Eq. (52) can be changed to a total 
derivative. Also, from the compressible Bernoulli equation, Eq. (52) be- 
comes 


Streamwise momentum equation (P = P(?) only) 


ff 


.. + [sm “ 
L a H - 


+ 3 




. 2 it- f sf 


- f 


^ d £n a H 


a in ? J/n? d fn § 


(55) 
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Energy equation 


•/* 


fH_ 



+ 


2 f 


3H„ 


3 in l 


- H, 


3 f 


T 5 in § 


(56) 


where q a * is the normalized diffusive heat flux away from the surface given 
by Eq. (50) and q* is the normalized radiant heat flux toward the surface 

q r * = q r /<X* (57) 

where a* is given by Eq, (44) . 


Elemental species equations 






(58) 


where is given by Eq* (48) . The transformation also yields the follow- 

ing relations between f w and P w v w (see Appendix B) : 


P 


v 

w w 


- a* 



+ 2? 


df 
w 

d§ 


(59) 


where a* is the normalizing parameter defined by Eq. (44) or, equivalently 


f 


w 



p v d§ 
w w 


p u a r 
*e e'e o 


(60) 


When certain groupings of parameters are constant so that the similarity 
assumption is valid, the terms on the right-hand side of the conservation 
equations (Eqs. (52) or (55), (56) and (58)) vanish, in which case the con- 

servation equations become ordinary differential equations. It should be 
emphasized that the equations as presented herein are equivalent to the 
corresponding boundary-layer equations presented in Section 2. That is, 
no similarity assumptions have been made in their development. 

Equations (44), (53), and (60) for a*, 3, and f , respectively, are 

w 

indeterminant at the stagnation point of a blunt body. Special forms for 
these equations valid at the stagnation point are shown in Appendix C to be 
given by 


a 


* 


o 



du 

< 

e ds 


3 


H 

o 


(44a) 


f w =-<Pw V w /a * ) o 


(60a) 
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where for Newtonian flow 


e o = i/u + i) 


and 


(53a) 



(2P/p)^/ R ef f 


with R eff an effective nose radius taking into account the shock shape. 
Alternatively, and (du e /ds ) q can be computed from curve fits of the in- 
viscid pressure distribution. 

In addition, in order to improve the accuracy of numerical integration 
procedures in the nose region, § and f can be computed by the following 
relations 


§ = 


2( 


Att/ 


f 

P u 
K e e 1 


lopl 
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Is/ 1 

s 1 
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d(s* H+8 ) 
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(h + 1) 


K + 1 
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p V 
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r o 

l s 1 

H- 

o 

- 


- 


d(s H+1 ) 


(31a) 


(60b) 


which take advantage of the fact that u q /s and r Q /s vary more nearly 
linearly in the stagnation region than do u g and r Q . Equations (31a) and 
(60b) are also derived in Appendix C. Of course, the°original Equations (31) 
and (60) are more applicable on the afterbody. 

The surface boundary conditions can assume numerous forms. The simplest 
of these are the requirement of zero slip at the surface which yields 


and assignment* of numerical values for p w v w (or f^ , h w (or , and K kw . 

In the event that p w v w is assigned, the can be calculated by use of W Eq. 

(60) . Alternatively, Eq. (61) can be utilized together with the assignment 
of wall mass diffusive fluxes, and h w (or Tj or with the assignment 


i«. P ^ ys i cally unrealistic in most cases to assign Kv when diffusion 
coefficients are unequal since the contribution to by preferential 
diffusion of the various elements to the surface is not known a priori. 
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of wall mass diffusive fluxes and the requirement that the surface material 
either be in equilibrium with the gas adjacent to the surface or satisfy 
surface reaction kinetic relations. (Surface chemistry considerations are 
discussed in Ref. 7.) 


Additional wall boundary conditions of interest admit the addition of 
chemically active species arising from the pyrolysis of an internally decom- 
posing material, surface combustion or phase change, and liquid-layer removal. 
In this case P w v w (and thus f^ by means of Eq. (60)), f^, H T ^ and 
are supplied through surface chemistry considerations, the zero slip condition 
(Eq. (61)), an energy balance, and elemental mass balances. 

The surface energy and elemental mass balances are supplied by transient 
internal conduction solutions such as those described in Ref. 20. The proce- 
dure for accomplishing this is discussed in Ref. 10. The resultant equation 
for the surface energy balance is given by 


m*h + m*h 
g g c c 


T. »; 

i 


( p v ) *h 
K w W w 


- q; 


w 


cond 


= 0 


(62) 


where is the mass flow rate per unit area and h^ the enthalpy of gas 

which enters the boundary layer without phase change at the surface (e.g., 
pyrolysis gases) , m g * is a normalized m^ given by m g * = m^/a* (typical) , 
m c is the mass removal rate per unit area and h c the enthalpy of surface 
material (e.g., char) removed by chemical reactions or phase change. ih is 
the mass removal rate per unit area and hg the enthalpy of that material which 
is removed in the condensed phase (e.g., by melting with subsequent liquid run- 
off or by mechanical spallation) . h^ is the enthalpy of the gas phase at the 
wall, q* is the normalized diffusive heat flux away from the wall (Eq. (50) 
evaluated W at the wall), q rw is the net radiative flux to the wall (including 
reradiation from the surface), and q conc j - ^ w (dT/dy) w is conduction into 
the surface material (with \ w the thermal conductivity of the surface mate- 
rial) . The elemental mass balances are given by 


m*K 
* g k 


m*X 
c c,. 



(P W V W ] *\ 



0 


(63) 


where the subscripts g, c, r and w and the asterisk have the same meaning 

and j* is the normalized diffusive net mass flux of elemental species k 
J k 

away fr$m the wall, given by Eq. (48) evaluated at the wall. 

To illustrate a simplified special case for this surface boundary condi- 
tion, consider the case of steady-state ablation of a homogeneous material 
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such as carbon, neglecting mechanical spallation and radiation absorption 
and emission with the exception of reradiation from the surface. In this 
case = 0, m^ = 0, = P w v w , q r = -a e W T W 4 # K^ c is unity for carbon and 

zero for nitroge^i, and q^ . = ni (ff - h ) where h = 0 if the carbon is 
initially at the base state of 298°K. Therefore Eqs . (62) and (63) become 

(p w v w ) * h w + q I + CTe w T w 4/a * = 0 < 64 ) 

w 


for the energy balance and 

<P w V*flfc - 1) + j* =0 

W W (65 

“V'J'V \ ' 0 

/ 

for elemental species balances for carbon and nitrogen, respectively.* In 
the absence of mechanical removal, equilibrium of the gas phase at the sur- 
face can be satisfied by considering any one of the equilibrium relations, 
for example 


(1/3) tn P c - K p (TJ = 0 

3 w eq 


( 66 ) 


where Pc is the partial pressure of C 3 at the wall and K (T ) is 
3 w P eq w 


the equilibrium constant for the reaction 


C (s) — ► ( 1/3) C - 


(67) 


Eqs. (64) through (66) together with Eq. (61) comprise the complete set 
of surface boundary conditions for this specific example of steady-state 
carbon ablation. 


*It is necessary to consider individual elemental species balances for only 
one less than the number of elements (see footnote at bottom of page 27) . 
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Returning to the general problem, edge conditions of interest admit the 
possibility of an entropy layer: 


edge 


h t 


= dr- 


edge 


edge 


edge 


^ (f ' ?) 


a H \ 


- edge 


df 




J 


edge 


edge 


actual 


= 0 


( 68 ) 


(69) 


V,„ = V 


edge 




edge 


eage 


actual 


= 0 


(70) 


where the subscript n e" refers to a reference condition, conveniently taken 
as the f = 0 streamline (see Appendix D) and the subscript "edge" refers 
to i1 e £g e chosen to be outside of the boundary layer but possibly in the 
entropy layer. When there is no entropy layer, the reference condition e 
can be considered as the edge condition, e = edge. In this case, Eqs. (68) 
simplify to: 


edge 


= f* = 
e 


H 


f " 

edge 


f M = 0 

e 


> (71) 


In the next section, the boundary-layer equations and boundary conditions 
presented in this section are cast into a form suitable for numerical solution 
by an integral matrix method. 
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SECTION 4 


THE LAMINAR BOUNDARY LAYER EQUATIONS 
IN INTEGRAL MATRIX FORM 


The solution of the transformed boundary layer equations presented in 
Section 3 utilizes an integral matrix method which has been developed spe- 
cifically for the solution of chemically reacting, nonsimilar, coupled bound- 
ary layers. In this procedure, the primary dependent variables f, H^, and 
K k and their derivatives with respect to r\ are related by Taylor series 
expansions such that f * , H^, and are represented by connected cubics 

with continuous first and second derivatives at the junction points (commonly 
called a spline fit) . Primarily for convenience, the conservation equations 
are integrated using a weighting function which is unity between adjacent nodal 
points and zero elsewhere. The linear Taylor series expansions together with 
linear boundary conditions form a very sparse matrix which has to be inverted 
only once for a given problem. The nonlinear boundary- layer equations and the 
nonlinear boundary conditions are then linearized, the errors being driven to 
zero using Newton-Raphson iteration. 

In this section, the Taylor series expansions are presented. The bound- 
ary -layer equations are integrated, and the integrals which appear are also 
expanded in Taylor series. The resulting equations are precisely those which 
have been programmed for solution on high-speed digital computers to repre- 
sent a coupled chemically reacting boundary layer such as surrounds an ablating 
heat shield during superorbital reentry. Special cases corresponding to a 
nonreacting (homogeneous) boundary layer and to an incompressible boundary 
layer are also discussed. The procedure utilized for solving the sets of lin- 
ear and nonlinear algebraic equations developed in this section are presented 
in Section 5. 

Consider the boundary layer in the region of a given station s as being 
divided into N-l strips connected by N nodal points. These nodal points are 
designated by r)^ where i = 1 at the wall and i = N at the edge of the 
boundary layer. 

Consider a function p(x) which with all its derivatives is continuous 
in the neighborhood of the point x = a. Then, for any value of x in this 
neighborhood, p(x) may be expressed in a Taylor series expansion as 


>(x) = p(a) + (x - a) + (x - a) : 


+ e ' - 3 -f- a) (x - a) 3 + P"^ (x - a) 4 + ... (72) 
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Considering the point a as 


r\^ and the 


as 


'i+l 


Pi +1 = p i + p i 6T1 + p'i + p'i" + p i' 
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( 73 ) 


where 


6ti 


'i+l 


(74) 


The p ± can be considered to be any of f^, f|, f£, f V' , H T , , K^, ! , R k , 

K k , or KT . To illustrate: 
i i 

f i+l = f i + f i 6t 1 + f i‘ ’ ^2~ + f i"' ^j“ + ••• (75) 

Since the highest derivatives of the dependent variables which appear in the 
boundary- layer equations are f| !J , Hj>\ and it is reasonable to truncate 

the series at the next highest derivative and to consider that derivative as 
being constant between T) i and 1# that is : 


f 1 1 1 1 

1 i+l 


fin . -phi 

i+l i 


6r] 


y III 

1 T i +1 


H * 1 - H ' 1 

T i +1 T i 


v II I 

i*k. 


i+l 


“ “k. 
i+l l 

6t] 


> (76) 


The following set of £3 + 2(1 + Kj](N - 1) linear equations are thus obtained 
where K is the number of elemental species minus one* and N is the number 
of nodal points across the boundary layer: 


f i+ i + f i + 


fpT, 


+ f i + f i ' + 


f... ini = 

r i+l 24 


(77) 


*For example, if 4 elements C, H, 0, N are under consideration, K = 3. 
Conservation of the remaining element is supplied through overall mass 
balance . 


27 


0 


(78) 


P i+ 1 + Pi + Pi 6T1 + p i + P i+1 ^6" 


- p i + l + p i + p i ¥ + P i+1 ll = 0 


(79) 


where in Eqs . (78) and (79) the p. represents f' represents H T , and 

^ 1 1 i 

represents each of the K . 

Thus, between each i and i+1 the f is represented as a quartic; 
the f ' , H t and are represented as cubics; and the f" , and are 

represented as quadratics; whereas the f"\ H^' and K£' are considered to 

vary linearly and the f "" , H^" and K£" are considered to be constants be- 
tween each pair of r)-stations. Of course, all of these functions with the 
exception of f ,,,, , 1 and K^" join continuously at the nodal points. 

Herein lies one of the major distinguishing features of integral methods in 
general and of the integral matrix method in particular. Conventional finite 
difference methods are generally based on a representation which is not too 
unlike the first term in a Taylor series expansion and thus yield discontinu- 
ous functions. Integral methods generally use smoother functions and hence 
can be made to yield comparable accuracy with far fewer nodal points. This 
is extremely important for a chemically reacting boundary layer since the 
state of the gas (the computation of which is not trivial) must be determined 
at each nodal point during each iteration. In the past (e.g., Ref . 3) inte- 
gral methods have employed high-order polynomials from the wall to the 
boundary-layer edge to obtain smooth profiles for f ' , H T and K^. The 
cubic spline functions employed herein are believed preferable as they are 
usually better behaved. 


The momentum, energy, and elemental species equations are integrated 
at constant § between and to yield; 

Momentum equation : 



»“h / T *’ - s / 
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Elemental species equation s : 


tK^ dr| 


- M* 


\dr) = 2 


r L 55l_.5e ,•« L 

/ r Un5 ^ 5 in 5 dT1 

i-l\ / 


This, in effect, is a square-wave weighting factor of unity between i-1 and 
i and zero elsewhere. This is equivalent to the step- function weighting fac- 
tor used by Pallone.^ As discussed in Appendix F, the primary advantage of 
this type of integration is algebraic simplicity, the complex terms in the 
energy and elemental species equations (the q*, q*, and j*) being divergence 
terms. The use of smoother weighting functions, such as those utilized by 

4 

Dorodnitsyn, would add considerable algebraic complexity and do not appear 
warranted on the basis of studies described in Appendix F which indicated 
that square-wave weighting functions and Dirac delta functions (i.e., a dif- 
ferential approach) yield comparable results as long as equivalent smooth 
functions are used to relate the primary variables to their derivatives. 


The integral of fp', where p is f, H^, or K^, can be expressed 


i r -i 1 i 

J fp' dr) - f p - f f'j 

*- i . i-i 


The integral of f’p can then be expressed by expanding in a Taylor series 
about i 
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f'p dr) = f: Pi ff - (f! P : + f: Pi ) f f- 


i- 1 


+ (f! p‘: + 2f'[ p ! + f"'p.) If 


- (f! P V' + p'| + 3f V' p^ + f!'"p.) j?- 

+ (4f: P V' + 6f'" p'' + 4f!-pl) 1 

- ( lOf p?" + lOf V" p" ) + 20f p!" 4IL + 

1 1 116. 1 r i 7 1 


(84) 


where, consistent with the truncation of the Taylor series employed earlier, 
all derivatives greater than fV M and p ” 1 have been dropped. Utilizing 


again Eqs . (76) to eliminate p V 1 , Eq. (84) becomes 


where 


f p dq 


i-1 


f! xp 1 + fl* xp 2 + fV’ xp 3 + f-; xp 4 


XP n 
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6 r) 


xp. - 6n 


XP„ 


( Pi - p£ I 3 * p;' + P iLi ) 
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(85) 


> ( 86 ) 


For the axial derivatives, logarithmic two- and three-point difference 
relations are utilized, namely 
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(87) 
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d In § 
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+ d 1 ( ) l _ 1 + d 2 ( ) 


je-2 




where ( refers to the previous streamwise station, and 

2 j 2 
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d l " " ,A 


l u l-X 


d 2 = 0 


( 88 ) 


for two-point difference and 


d 0 =2 


-e A i-i + i A je-2 

jC A jC-2 


- 2 


l L l- 2 




a _ n 
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i A X-l 


f i-2 X-l u i-2 


> (89; 


for three-point difference where typically 


X A / - 1 “ in “ 411 5 i_l = in 


(90) 


The three-point difference relation is utilized unless a similar solution is 
desired (in which case d Q = d^ = = 0) or unless the point in question is 

the first point after either (1) a similar solution or 2) a discontinuity 
(e.g., where the body changes shape abruptly, or where mass injection is sud- 
denly terminated) . 

Similar approaches have been utilized previously in finite difference 
procedures, for example, by Smith et al and by Leigh, Integral methods, 
on the other hand, have generally integrated the boundary-layer equations in 
the streamwise direction (e.g., Pallone^) . The present approach is considered 
preferable since the streamwise derivative terms are generally small, and 
hence are conveniently treated as forcing functions. This avoids the diffi- 
culty sometimes experienced during integration when the streamwise step size 
is too small (see Ref. 6) . 

Applying Eq. (87) to the streamwise derivative terms 
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( 91 ) 


where again this relation applies for the three cases where p is equal to 
f,H T and 1^. Utilizing Eq. (83) ‘yields 
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(92) 


Noting that each of the five integrals in Eq. (92) is of the same form as the 
integral expanded in Eq. (84) by means of Taylor series 
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(equation continued on next page) 
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(equation continued from previous page) 


+ 2d o (f i XP 1 + f i' * p 2 + xp 3 + H-1 * p 4 > 

+ (f[ ZP 1 + f*' ZP 2 + fj" ZP 3 + f ZP 4 ) 

+ (p. ZM X + p| ZM 2 + p» zm 3 + pM^ m j 
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and ZM^, ZM 2 , ZM^, and ZM^ are equal to ZP^, ZP 2< ZP^, and ZP^, respec- 
tively, for the special case that p = f 1 . 
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Following the same procedure 
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-J p. = f! 

*1 l 


where use has been made of the fact that the coordinate stretching parameter 

a„ is a function of § only. 

H 

Finally, it is necessary to evaluate the integral of the density ratio 
which appears in the momentum equation and the integral of the elemental 
source term which appears in the elemental conservation equation. Approximat- 
ing P e /P as a cubic between i and i-1, an exact integration of the result- 
ing approximate integral yields 
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integral of 0^. 
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(97) 


(98) 


These approximations are not quite as good as the approximations for f ' , H fJ 
and since continuity of derivatives is not guaranteed at the nodal points. 

Utilizing Eqs . (83) through (98) the boundary -layer equations (Eqs . (80) 

through (82)) become 
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(equation continued on next page) 
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(equation continued from previous page) 
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Energy ; 
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( 100 ) 


and q* are given by Eqs . (50) and (57), respectively. 
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Elemental species ; 
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( 101 ) 


where j£ is given by Eq. (48) . 

The boundary- layer equations (Eqs. (99) through (101)) are applicable to 
the problem of the nonsimilar chemically-reacting boundary layer with unequal 
diffusion coefficients, thermal diffusion, entropy layer, radiation absorp- 
tion and emission, and rate-controlled reactions, coupled point by point with 
a charring ablator solution. In the absence of thermal diffusion c t is set 
equal to zero. When the diffusion coefficients are equal, = K^, = 1, 

4 ^ = 1/^3 = Wit h = h, and = C . When the boundary layer does not react chem- 
ically, the elemental species equation is inconsequential. Finally, if the 
boundary layer is incompressible or if Crocco relations are utilized only the 
momentum equation is needed. 

Before discussing the procedure for solving the equations developed in 
this section, it is appropriate to discuss briefly the thermodynamic and 
transport properties employed in the solution procedure. These subjects are 
treated in considerably more detail in Refs. 7 and 9, respectively. 

The state at each node is determined with a general purpose chemical 
equilibrium subprogram of much the same form as those described in Refs. 21 
and 22. State derivatives are determined by the same routine. Enthalpy and 
specific heat values are obtained through accurate curve fits of JANAF or 
other reliable thermochemical data. 

In the present formulation, transport properties are determined as fol- 
lows. The and are calculated using the approximate Eqs. (19) 

and (26), respectively. The D is given by 
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D 


2.628 x 10" 3 


T( T/^ref) h 


Pa; 


ref ref 


(1,1)* 


(cm 8 /sec) 


( 102 ) 


with T in °K, P in atmospheres, and a in A. The subscript "ref n refers 
to a reference species (often 0^$ but conceivably fictional) . D is thus the 
self-diffusion coefficient of that species. The are determined by a 

least-squares correlation of 3^ - 

The viscosity of the mixture is obtained from an approximate relation of 
the Sutherland-Wassil jewa type 13 




where j-u is the viscosity of species i given by 


6A ii p ^ u 


(103) 


(104) 


and 



1 + 1.385 


RTH i 

PxpiT 



(105) 


with the 1.385 an empirical constant suggested by Buddenberg and Wilke. 

The mixture thermal conductivity, is obtained as the sum of a mona- 
tomic thermal conductivity and a contribution from internal degrees of freedom: 

X = X mono + x int (106) 


The \ 


mono 


is determined by a relation similar to that for \i 


^mono _ ^mono (p 


(107) 


where \ mono ± s the monatomic thermal conductivity of species i 

i 


. mono 
A i 


4 ^ 


(108) 
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and 


(? 


1 + 1.065 x 1.385 


RTn i 

PxT^T 



(109) 


with the 1.065 being suggested by Mason and Saxena.^ The is computed 

25 

using 


x int 


l 


x i (X i - 



. mono* 

ii 1 



where 


X . 
i 


^mono 
* i 


\ 

p ^ii 7n 




SECTION 5 

SOLUTION OF THE BOUNDARY-LAYER EQUATIONS 
IN INTEGRAL MATRIX FORM 


( 110 ) 


( 111 ) 


The solution of the boundary-layer equations presented in Section 4 to- 
gether with the boundary conditions such as those presented in Section 3 is 
accomplished by Newton-Raphson iteration. In this section these equations 
are put into a form suitable for solution by this procedure. The resulting 
equations are then written in matrix form, and a method is presented for their 
solution suitable for coupling with an internal conduction solution. The pro- 
cedure attempts to minimize computational time and computer storage require- 
ments . 

In order to illustrate the Newton-Raphson method consider two simultane- 
ous nonlinear algebraic equations 

F (x, y) = 0 G(x,y) » 0 (112) 

the solution for which is given by x = x, y = y. D 3 fine x and y as 

th m m 

the values of x and y for the in iteration. The desired solution 

f(x#y) can be expressed in a Taylor series expansion 
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0 = F(x,y) 


+ < x " x ™> 




m' "53T 


m' ■‘nr 


dF(x ,v ) 

+ ly _ v ) + 

ly ”m' dy + . . . 


) ( 113 ) 


0 - G (x ,y) - 


3G(x ,y ) 

G( x ro ,y m ) + (X - Xj 


m Ax 


dG(x ,y ) 

+ (y - y ) ^ — m + . . . 

Ay 


The Newton-Raphson method consists of replacing (x, y) by (x , , y , ) on the 

m+ i nw - 1 

right-hand- side of these expressions, neglecting nonlinear terms in x " x m 

and y - y . This yields the set of recurrence formulas 
m+ i m 


Ax 


3F(x ,yj 


m 3x 


+ Ay 


3F(x,yJ 


m dy 


- F(x m>y m ) 


AG (x ,y ) AG (x ,y ) 

Ax_ — 1 SIISL + A y - m m 


m Ax 


Y m Ay 


- G(x,y) 
m'^m 


l (114) 


J 


where 


Ax = 


m 


x , - x 
m+ 1 m 




3 Y, 


m+ 1 


- Y, 


m 


(115) 


The Ax and Ay are the corrections to be added to x and y . respec- 
m m mm 

tively, to yield the values of the dependent variables for the m+i^* 1 itera- 
tion. Here F(x ,y ) and G(x ,y ) are the values of the original functions 
m m mm 

F ( x , y) and G(x, y) evaluated for x m x^ and y ■ y^. A9 the corrections 

approach zero, the F(x ,y ) and G(x .yj thus approach zero. Hence, it is 

appropriate to look upon these as errors associated with the original Equations 
(112) . It is apparent that this procedure can be extended to an arbitrary num- 
ber of functions and a corresponding number of primary variables. 

For the purpose of the present problem, it has been found to be most con- 
venient to consider the primary variables as the f^, f^, fV , fV' , H T * h t . ' 
HrJA, K k ., K k ., K k , at each nodal station i, plus the a H - This amounts to a 

total of (7 + 3K) N + 1 unknowns where N is the number of nodal points across 

the boundary layer and K is one less than the number of elements present in 
the boundary layer. The corresponding number of equations is provided as follows 
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Taylor series expansions 
Boundary layer equations 
Boundary conditions 

The a u constraint 

rl 

Total 


Equation numbers 

(77)- (79) 

(99) -(101) 

(61) - (63) and 
(68) -(70) or equiv. 

(34) 


Number of equations 
(5 + 2K) (N - 1) 

(2 + K) (N - 1) 

( 7 + 3K) 

1 

(7 + 3K) N + 1 


It is thus necessary to specify the corrections of the other variables (such 
as density, temperature, etc.) in terms of corrections of the primary variables. 
The procedure for accomplishing this will be described later. 

The Taylor series expansions are linear with respect to the primary vari- 
ables as are several of the boundary conditions. The boundary layer equations 

and the remainder of the boundary conditions are nonlinear. The a constraint 

H 

is linear but it must be considered together with the nonlinear equations in 
order to avoid a singular matrix. The recurrence formulas representing the 
linear equations will be presented first, after which recurrence formulas 
appropriate to the nonlinear equations will be developed. 

Partial differentiation of the Taylor series expansions with respect to 
the primary dependent variables in accordance with Equations (114) yields for 
the m th iteration 


(-1) Af i+1 + (l)Af. + (&n)Af[+ 

(116) 

(117) 

(118) 

and 
the 

numbers in parentheses represent the partial derivatives of the Taylor series 
expressions (Equations (77) through (79)) with respect to the primary variables; 
and the ERRORS are obtained by evaluating the left-hand-sides of the appropriate 
Equations (77) through (79) for the values of the variables obtained during the 
m t ^ 1 iteration. 

Similarly, the recurrence formulas for the linear boundary conditions 
( Eqs . (61), (69), and (70)) are given by 


+ (^-)Af(" + (|f = - ERROR 

(-!)APi + l + (D^ + (6r|)Ap^ + (_^Dl)Ap^ + - _ ERROR 

(-D4p| +1 + (l)Ap[ + (y^ApV + (y^Ap^ - - ERROR 

where as before p^^ represents f^H^ and K k . . Here Af i+1 > Af^ A f^, 
so on represent the respective corrections for f^ + ^, f^> and so on > 
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Af 1 = 

w 


AH,, 


- ERROR = - 


edge 


<f') m 
w m 


= - ERROR = 


■h 


- H m 
T 

edge edge 


actual - 1 m 


( 119 ) 

(120) 


AH' 


= - ERROR = - 


edge 


(H' ) 

v T - ' m 

edge 


(121) 


AK k 


= - ERROR = - 


edge 


\ \ 

edge edge 


actual - 1 m 


( 122 ) 


AK£ = - ERROR = - (K/ ) 

edge edge m 


( 123 ) 


The recurrence formulas for the nonlinear boundary-layer equations are 
given by: 

Momentum : 


C f‘ 


Af n AC 

f" + C 


Aa 

sf ] * | (1 * d 0>f + Vi-1 + Vi 


- 2 ] 


f' (1 + d Q ) Af 




-I i-1 


H p ± » 2 


3 P ± 


APf - ^ Ap<_ 


"i 

Pi-1 


_ 6tl fUl 
3 Pi-1 


Ap 


i-i + t 1 Ap i-i 


+ ea H 6 ^ 77 


1 + 


Pi-1 


in 

6 







,d,a TT + d~a TT 

f Pi 

pi p i-i i 


Aa H - 

1 + 3 + “| 

1 H £-l 2 H f -2 

IpI' 

pi-i p i-i j 


i I] 


f! AXP. + f'.' AXP_ + f!’ 1 AXP, + f. M, ,AXP. + XP,Af! + XP,Af!' 
l 1 l 2 l J i-l 4 li ^ i 


+ XP 3 Af i M ' + XP 4 Af 




d i a H + d 2 tt H 

1 H it-1 2 jt-2 


' £ i ' 


f :xp. 

i i 


(equation continued on next page) 
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(equation continued from previous page) 


+ tr XP 2 + f«" XP 3 + f ! 


Aa H - 2 

Pi = q 


+ ZP 2 Af[’ 


+ ZP 3 Af|" + ZP 4 Afj 


= - ERROR 


Pi - q 


where the ERROR is given by the left-hand-side of Eq. (99) evaluated for the 
th 

m iteration. 


- Aq* + Aq r + (1 + d Q ) f + djf^ + d 2 q_ 2 AH T + 1^(1 + d Q ) Af 


- (1 + 2d 0 ) £fj AXP^ + fV AXP 2 + f»* AXP 3 + f!^ AXP 4 + XP^q 


+ XP 2 Af V + XP 3 Af V 1 + XP 4 Af] 


Pi = H T, 


’ ZP l Af i 


+ ZP 2 AfV + ZP 3 Af|" 


+ ZP 4 Af!'J 1 J p _ H - L ZP 1^ H T . + ZP 2 AH Ti + ZP 3 AH t. + ZP 4 AH t. 1 J 


ERROR (125) 


where the ERROR is given by the left-hand-side of Eq. (100) evaluated for the 
m iteration and Aq* is given by 


Aq* = 


C f'f'u 3 


AC 

Ml 

. Af" 

3 ^h| 

CC T' / 

i p 

c H 

r ~ 

f " 

a H ' 

a H Pr 1 


c -p T' 
C p 


T-J - - —■ + — K'- c + — - — T' + c RTq' + (K - h 

a H Pr / a„Sc P ^2 I 3 


(equation continued on next page) 
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(equation continued from previous page) 


+ c t RTn 3 )n 4 



AC 


_ ASc 

c 

1 

% 

i 1 

C +- R ^ 

r a. 

. 

C 

_ 

a H 

Sc _ 

a H Sc 

P ^2 j 


AT' 


-T'AC + 

P (M 


t RT ’ / AT A ^3 1 

- 3 ^1^ + C t R ^3 \f + HI 2 + (h “ h + c t R^ 3 ) 

M 1 M 2 ) 2 ' M 3 / 


+ M 4 Ah - Ah + c t RTp 3 


AH 


3 + AT 


t 


Elemental species; 


Aj* + 


(i + d Q ) f + d^^ + d 2fji _ 2 


A \ + R k (1 + d 0 5 Af 


J i-1 


+ a 


in 

H 2 


h \ + A 0 k. _ x - I 1 [ A 0 k ± ~ A 0 k ._ 1 


+ in 
2 


\ + V 


*i. - *i. ,, 
1 1 - 1 / 


6n 

6 


Aa T 


+ f PJ 1 AXP 4 + XP 1 Af[ + XP 2 Af[" + XP 3 Af[" + XP 4 Af 


(1 + 2d Q ) £ f ^AXP^ + fV AXP 2 + f!"AXP 


pi - 


- r-i A - + ZP 2 Af i + ZP 3 Af i' ' + ZP 4 Af i-l] ~ " [ ZP l A *k. 

L J Pi = K k. 1 


!P A 1 + zp 3 ‘% * . f , ’ 


ERROR 


where the ERROR is given by the left-hand-side of Eq. (101) evaluated 
the m t * 1 iteration and Aj£ is given by 


^4 


(126) 


1 


3 


(127) 

for 
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Aj* 


_C 

d H S^ 


+ 


<\ - «k^4 


AC 

C 


-M\ 

a H Sc I 


+ AZ^ + (Z k - \)A^ + M^(A\ - AKj^) 


(128) 


Equations (124), (125), and (127) are reduced to linear equations in terms 
of the corrections on the primary variables (Af^, Af|, and so on) by noting 
that the variables C, p, C T, Pr, Sc, h, C p , V 3 * 4 4 , q r , Z k# and 

evaluated at any point in the boundary layer can be considered as func- 
tions of static enthalpy, static pressure, and elemental composition. With 
the pressure assumed constant across the boundary layer, it follows that all 
of the corrections on unprimed variables with the exception of Aq r can be 
expressed as 


M> i ■ Z 


»( h 


an 


AK kk. 


kk 


kk. 


3 ( ) 

ah“ 


Ah . 


(129) 


where from Eqs . (9) and (36) 


h. 


l 



2 «■ 


(130) 


so that 


Ah. 



f'. 3 

1 

2 

H 




(131) 


The Aq r is more complicated in that it depends upon the AK kk _. and Ahj 
at all nodal points j . 

The ^-derivatives of these variables (i.e., the primed quantities) can 
likewise be expressed in terms of corrections on the primary variables as 
follows 


A( )[ = 




kk 


). 

~ ~ AK kkk + "13 — Ah i 


+ 



5 2 ( ). 3 s ( ) • \ 

AK. t + Ah . 

dh i 9K kk. ^ h 



kk 


SK kk. 


AK k ki + Mr 1 Ah i 


(132) 
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where 


so that 


h. 

i 


H ' - 
T i 


i 2 f f ' 
e i i 


(133) 


u 2 f' f'.' 

. e r l 

Ah i = ah t. “ “ 7 

a H 


Af’ Af 


Aa T 


fT 

i 


tv 

i 


(134) 


Use is also made of the following which are obtained by differentiating 
Equations (86) : 

AXP X = 6ri | Ap. - f 1 Ap^ + Ap^' + A Pi_l) 

AXP 2 - - 6n“ ( - f- Apl + ApV + # A P i.J 

1 - SP API - ^ AP-' + Mf **u) 


(135) 


AXP 


AXP = 6r) 

4 


Ap i 

8 


Ap i 6r| . 5 6r] 3 ,, 6ri 2 * • » 

24” “ 30 Ap i + 504 Ap i + 252 Ap i-1 


The AZP^ = AZP 2 = AZP 3 = AZP 4 = 0 since ZP^, ZP 2 , ZP 3 , and ZP 4 can be com- 
puted before the iteration commences. 

Substituting Eqs . (126), (128), (129), (131), (132), (134), and (135) 

into Eqs. (124), (125), and (127), and collecting terms yields, neglecting 

the off-diagonal terms* arising from the Aq , the following recurrence forjnu- 

“* 4-"h * 

las for the momentum, energy, and K elemental species equations between 
nodal points i and i-ls 


c l Af i + c 2 Af i-l + c 3 Af i + c 4 Af i-l 


+ c 5 Af!' + c 6 Af!'. 1 + c ? Af!- + CgAfV^ 


+ c 9 AH Ti + c loA H T i 1 + C ll AH T i + C 12 AH T i _ 1 


(equation continued on following page) 


*As mentioned previously, Aq r ^ and Aq r ^_^ produce corrections AK^ j , 
AHipj , Afj and Aa H ^ at all boundary-layer nodal points j. 
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+ c 13 AH^ + c 14 AH T i _ i +c l5 Aa H i + C 16 ACt H i _ 1 

+ Y ( c kkl AR kk. + c kk2 AK kk. . + c kk3 AR kk. 
kk ' 1 1-1 1 

+ c kk4 Al? kk i _ 1 + C kk5 A ^k. *■ c kk6 A ^kk i _ 1 ) = " ERR0R < 136 > 


where the coefficients for the corrections c, , c 0 , . .. differ in form for 

the momentum, the energy, and the K elemental species equations and in gen- 
eral yield different numerical results for each elemental species kk and 
for each nodal position in the boundary layer i. The coefficients for the 
corrections are listed in Tables I through III for the momentum equation, the 
energy equation, and the K^ 1 elemental species equation, respectively. 

In order to complete the set of equations, it is necessary to develop 

the recurrence formulas for the a„ constraint and for the nonlinear boundary 

Jti 

conditions. The a„ constraint (Eq. (34)) yields 
rl 

Af ' - cAf* = - ERROR = - (f * - cf * ) (137) 

^e ^c ^e m 


The boundary conditions at the boundary-layer edge are nonlinear when an 
entropy layer is taken into account. In this case, Eqs. (68) are applicable, 
resulting in the following recurrence formulas; 


-Af edge + [u # (f ' §) ] edge Aa H 


+ a T . 


df 


Af 


edge 


- ERROR 


edge 




edge * »H [s“ (t ' 5 ^ 


edge 


J m 


(138) 
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TABLE I 

CORRECTION COEFFICIENTS FOR MOMENTUM EQUATION 


Af. : 

[ (1 * V f, ] 1 

4f i-l '-[» * a 0> 4 '] 1 _ 1 

Af ^ : 

M i * B i 


Af" : 

i 

■f-l * b 2 

- H-l i 2 

Af!' .. : - [SJ\ + c. 

1-1 L a H-l i-1 2 

Af!” : 

B 3 

Ar lll _ 

Af i-1 : C 3 

: 

T i 

[C w] 4 * <■>,>!_ . fh 

“*1-1 : ‘ fe i.,. + (D 4>|_ 

dp 

ah' : 

T i 

°3 llfeli 

aha : D„ (|£-1. . 

^ 3 ' on * JL-1 


th 


Ah' ' : 0 

T i 


AH' 


T i-1 


0 





* 'Vi. 



f 11 ac 

a H bK 


kk 




-j i-1 


bp 


b_ 

**kk 





Aa 


H i 





i-1 



0 


Aa 




i-1 


+ C 


4 


f f'f" !§ + d + V f + *1^-1 + d 2 f *_2 




■ - 2 fo, XP. + ZP,"| 

L 1 1 IJ p. - f'. a » 


f i °3 


l£ 

dh 


+ f i . a 

■§p '5h_ 
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TABLE I (continued) 


- 2 [°1 “2 * ZP 2 ] p . - f. - £ ‘‘M ^l 1 


Pi - h «; 


2 D ] _ XP 3 + ZP 3 


Pi = 


d,a + d.a 

1 H X-1 2 H jt- 2 


f [ XP X + f V XP 2 + f V ' XP 3 + f XP 4 


* -7 2t i n i - f i < d 2> s , » 


bp dh 


+ 3 a H 6r| 


p e 

-g. i . 

P i . M 

[ p i 

p i p l-i[ 

p i L 

r t r 

p i-i 6 1 

l p i 

p i-l p i-l /. 


fV , dJI^L , + f- (D d ); 


1 u 5 3h i-1 T x i-l ' u 4'b_ m b_ 
' 1 bp bh 


I® f < D l£ 

3 l-l U 5 3h i-1 


2 Dl XP 4 + ZP 4 


Pi - f i 


~f f i-l 2f i-i° 5 (fh)i-l + f l-l (n 4 ) S_ . 3^ 
a H L ' 3p 3h 


1 + 3 + d Q - 


d,a + d,a 

1 H X-1 H X-2 
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TABLE I (concluded) 


8 a 2 -#• 
P a H 2 p| 


/ 

n- \ 

/ \ 


i+ in fi) 

Ifi. 

6ti 

ir 3 °i/ 

Mi 

6 


E 

kkk 


K£ kk I ]+ 1»* 

*P 3K kkk 


lie. 


dp3h 


D, = 


e a 3 ini 

P a H 12 p| 


3 a: 


6r| H e 

H 2 Pl _ x 


_ in p l-i 


3 p. 


i-l 


ni 1 +± & r z ^ 

/ 1-i L kkk 


He. 


3 P 9 \kk 


+ h 'l life 


i-l 


°5 = - e a H TT ^f_ x 
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TABLE II 

CORRECTION COEFFICIENTS FOR ENERGY EQUATION 


Af i : 

(1 + 

d o lH r] , 

Af i-1 5 

-[»* 

*0> H t] 

Af'. : 
1 

[a 1 

L 15 J 

+ B 

i 15 

Af U : 

- M 

i-1 

Af'.' : 

1 

r A n 

L 16. 

+ B 

i 16 

* 

- [*J 

i-1 


, _ 1 1 1 „ 

Af. : B_ 

i 17 


*Ti : [ A ie] t * b 18 

i T. : ^° 14 ^ 3 9 

1 [ a£ ah J i 


'i-i ‘ -Wi-i 


- < D 14 ) 9 


- 

- 

31 ?J i-1 


AH T. : B 20 


Ai W <^ 3 >a _ a_ 
sic 




A *kk. " {D i3h m a 

* ^ , 


iK kk.' ln L4»S .9 

SP **kk , 


A *kk. " (D L4 ) 9 „ a 

Sp >** , 


0 

Sa H t - : [*21 ] i 


0 


■i-i ; -Wi-i 
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TABLE II (continued) 
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TABLE II 


(continued) 


*17 


a ♦ [xp 4 ] ^ _ * [» 4 ] 


p ■ - H m 
1 T . 


"20 


11 + 2a °’ H P . . £ ; * [ 2P <] 


1 1 


Pi = f ; 


j n 


1 2 


I C 3 R \ 

h’ - C + — — 

\ P 


c_ 

a,. 


«*H PC “ 


t' + c t R T M 3 + (h - h + c t R T P 3 )M 4 

+ < j £ 


_ „** I 3 £ . 1_ 4. _C 

d 13 “ (D 12 “ q r> C 5p + a* 5p + a H Pr 


!fs . iit 

3p Pr 3p 


a H Sc 


(P 4 P 3 + M 3 )c t R If + U 4 Ip + M 4 C t R T Tp" ~ T ' Tj? 


Sc 


3S? 1 c t RT ' 


C M . 

Sp ( Ml u 2 )» 5p 


_ a H Sc 4 . 


8 p = 9 h onl * 


C C 


£ _ 


, c + 


G t R 


a H Pr a H s5 \ P ^2/ 


L 


kkk la P 9K kkk 


9 S T 


+ h' ll3L 

| Spdh 


t 

a H Sc 


(c. R T) 


E 

kkk \ dpdK kkkl 


+ h 


s 2 u 3 

\ 5p§h 


a H Sc 


(h - h + c fc R T n ) 


I ^kk 




kkk 


dpdK 


+ h 


kkk i 


. I 

dp3h 
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TABLE II (Concluded) 


_C ^ ST + C 

3 14 = a H Pr p •§? a 

H 


Sm, du 

c t R T TT + (h - h + c t R T u 3 ) ^ 


Is + °^] 

ST 

4_ 

C 

(p 

3p 

T 

-“h 55 - 


Tp Sh 


only 
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TABLE III 



: 

i 

w 

dp dh 

“Vi : 

-H 

j-i, a , a_ 

5p 5h 


- a H • 

firf. f a *k] 
12 [lirji 


+ a H 

±£. [!!k| 

12 Uhl i-l 

Ah t. : 

0 


: 

0 


A *W 

1 

[v 

i + B “ 

s ^i-i ; 

- |V 

i-l + Cn 

A ^k. : 

M 

L* ^ 

si ^i.i : 

- [* 12 ] 

+ a 

i-l H 12 

A ^k. J 

=13 


~n 

‘^i-i : 

C 13 


Aa H. : 

M 

♦ =14 

i 

la »i-i : 

- [v 

. + C 14 

i-l 14 
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TABLE III (continued) 


A 7 “ 


— ( f " D c + f D ) 


'6 T " 7'd _ d 


'll' Hi 9 + 

■ 


(1 + d 0 )f + d^J, + d 2 f t 


-2 " -V ^1 

a H SC -I 


k ■ kk only 


BK 


kk 


A 12 L D 6j a a 


Bp 


BK 


kk 


*14 - £ *’ < 2f " ft * ^ - *>"*] 

H op Bh h 


- - ^ZP X + (1 + 2d 0 ) XPjJ 


Pi ■ ^ 




f . _ ±a f..| [i^. 

f i 6 i I dh J 


i - T 1 f i |D S> i. 9 9 

Jp = 9h* 


- ZP 2 + (1 + 2d 0 ) XP 2 


u e 6n 3 f . I 5 

n - K. a H “TTT f i 3h“|i 
P i % a H 1 


B 9 “ - [ ZP 3 + (1 + 2d o> xp a] 


Pi “ Kk, 


B 10 = a H ¥ 


8 *k 

aTj 


i " T 1 (D 8 J i. 3 _ 3 

3p ' 3h- 
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TABLE III (continued) 


= - ZP-l + (1 + 2d Q ) XP 1 


Pi = f|, k = kk only 


+ a H tF ^ - T 1 ( V i, d_ _ 1 

l “kk/ i 3p = ^ 


B 12 ■ - [ ZP 2 + 11 * 2d 0> =*2] ■ a H if (-rH 

L J Pi * f i- * - ** W\ k /i 


= - ZP 3 + (1 + 2d 0 ) XP 3 


= f|, k = kk only 


!®M f . Nf. _ la f ..l _ m f . ( D ) 

n 3 2 T i I r i 3 r i J \ dh / i 6 1 (D 8 J i, 3 


* f [\ * Vi ■ K ■ \.J f] 


a Jjain f . + An f « _Js + in f. > 

H a 2 | i-l + 6 i-ljldh / i-1 + 6 f i-l (D 8 i-1, 9 d_ 

H L. dp dh 


a H a 12 f i-l l dh / i-1 
a H 




ZP 4 + (1 + 2d 0 ) XP 4 


Pi = Kk. 


C 10 ~ a H 


la Tf!^) + la (D ) 

2 dh | i-i + 6 i-i, a_ a_ 

L dp dh-l 


Cn= 
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TABLE III (concluded) 


'13 


ZP 4 + (1 + 2d 0 ) XP 4 ] 


P i = f|, k = kk only 


c _ _e in f , 

C 14 a 2 f i-l 

tt H 


fi + M fn 
r i-l + 3 i-l 


* + in f i ( D ) 1 

9h j i-1 + 6 £ i-l D 8 i-1, 9_ _ dl 

Sn ~ ShJ 


3p dh- 


°6 = "^= 
a H Sc 


Ifk + _ V ) ^4 

dp ^ 3p 




a H Sc 


. *** + 
> — + 
dp 


ft 


dSc 

bp 


Z! *kkk 


kkk 


* Z k _l +h . 

+ h 1 bpbh 


3 P 9K kkk 


+ - *k> 


E 

kkk 


r (!% ) + h >f d % 

K kkk 1 ) + h Up§h 


\ d P di W 


D„ = 


E ^kkk 


9 2 0, 


kkk 


i9 P 9i W 


u 


dpdh , 
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Af - + 

edge 


fe (f - 5) ] 
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(139) 


In the absence of an entropy layer, Eqs. (71) yield the following recurrence 
formulas : 

Afg - Aa H = - ERROR = - (f^ - a H ) m (140) 

Af" = - ERROR = - f" (141) 

© s 

Although Eqs. (140) and (141) are linear, these boundary conditions are in- 
cluded in the nonlinear set because of the possibility of an entropy layer. 

As pointed out in Section 3, the wall boundary conditions can assume a 
variety of forms. When the boundary layer is coupled to a transient charring 
ablation solution, there are a total of K + 2 nonlinear wall boundary con- 
ditions supplying, in effect, P w v w (and thus f^ by means of Eq. (60)), 

Hfp w and K Kk w - Therefore, although some of these boundary conditions 
become linear for simpler problems, they also must be considered as nonlinear 
for the general solution procedure. Solution for the wall boundary conditions 
for the fully coupled problem are discussed in Ref. 10. The wall conditions 
for two simpler problems are described herein. 

In the case that the surface conditions are assigned a priori, these 
boundary conditions become particularly simple. For example, if P w v w ' 

% w and K kw are assigned, the following recurrence formulas result: 
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d (p v ) * 
w w' 

df 


Af w = - ERROR = - | (P w v w )* - (p w v w ) 

[ h t - h t 

I W W 


ERROR = -I 


actual 


actual j 


m 


1 

** 

1 


I w w | 

actual 


where from Eqs. (59) and (87) it follows that 

a( Pw v w)* 


3f. 


(i + d Q ) 


(142) 


(143) 


(144) 


(145) 


To illustrate a more complicated situation, consider the case of steady- 
state ablation of carbon. The appropriate wall boundary conditions (given by 
Egs. (64) through (66) ) yield the following recurrence formulas: 


(p w v vP * AI Y. + V. y r w - Af w + Aq S +_ st Lt w 3 AT w " " ERR0R (146) 


w w w 


w 


d (p v ) * 

(p w v w } * AR C + (R c " 1} af W W Af w + A ^c = - ERR0R 
WWW w 


> (147) 


(p w v w } * AK N + ^ df 


— P w V w?I Af w + Aj* = 


ERROR 


w w w 


3 A(Xn p c > " W 


9K_ (TJ 


_eg_ 


AT = - ERROR 

w 


(148) 


w 


where the ERRORS are given by the left-hand-sides of Eqs. (64) through (66), 
respectively, evaluated for the iteration, the is given 

by Eq. (145), the Aq * ' and are given by Eqs. (126) and (128), respec- 

tively, evaluated at the wall, and the AT w and A (in Pq ) are reduced to 
the primary variables by use of Eq. (129). w 

In many problems it has not been necessary to consider all of the terms 
in all of the coefficients in the recurrence formulas. The retention of the 
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major terms has been seen generally to improve convergence and thereby reduce 
the number of iterations. However, overall computational time and storage 
requirements can be improved by dropping the lesser important terms. At pres- 
ent, the terms involving second derivatives and some of the lesser important 
terms involving first derivatives are excluded. Of course, the dismissal of 
terms from the recurrence formulas does not affect the accuracy of the final 
result as long as the ERRORS are evaluated precisely? rather, the errors are 
driven to zero along different paths. 

The coefficients for the recurrence formulas for the Taylor series ex- 
pansions (Eqs. (116) through (118)), the linear boundary conditions (Eqs. 

(119) through (123)), the constraint (Eq. (137)), the nonlinear edge 

boundary conditions (Eqs. (138) and (139)), and the boundary-layer equations 
(Eqs. (136)) evaluated for the III th iteration form a non-square matrix Qa]] 
with I = (7 + 3K) (N - 1) + 6 + 2K rows (the number of equations, excluding 
the nonlinear surface boundary conditions) and J = (7 + 3K) N + 1 columns 
(the number of correction variables) . This matrix equation is given by 

[ A ] j^AvJ = ~ [ E ] ( 149 ) 

IxJ J I 

where AV • represents the correction on the j th primary variable (Af,, 
j th 

Af , ..., Af^, etc.) and E i represents the error associated with the i 

equation . 

This matrix can be reduced such that all corrections are expressible in 
terms of Af , AHm and the AK* . This approach makes it convenient to 
treat varied and complex surface boundary conditions. Any consistent set of 
surface boundary conditions can be added as an option with a minimum of pro- 
gram modification. The primary reason for this approach, however, was to im- 
prove the numerical stability if the boundary layer were to be iteratively 
coupled to a transient charring ablation solution procedure, by relegating the 
highly nonlinear surface boundary conditions to a subsidiary iteration with 
the charring ablation solution. The influence of the boundary layer would be 
contained in the reduced coefficients of the Af w , AH^ an< 3 ^ K k w somewhat 
analogous to convective transfer coefficients. 

The matrix [a] is quite sparse (i.e., it contains many zeros) in an 
orderly way. Substantial savings in computation time and storage allocations 
can be realized if full advantage is taken of this ordered sparseness. This 
is extremely important since the solution of a boundary layer with several 


60 


elemental species would otherwise be very costly. For this reason, the matrix 
solution procedure will be discussed in some detail. 

The first step in the matrix solution is to divide the equations into 
linear (symbol L) and nonlinear (symbol NL) sets, namely 


AL 

! BL 
1 


AVL 


EL 

ANL 

l BNL 


AVNL 


ENL 

_ 

1 




_ _ 


(150) 


I x J 


where for convenience the variables are also classified as "linear" and "non- 
linear" . The distribution into linear and nonlinear variables is somewhat 
arbitrary, but care must be taken that the square matrix Qal] not be singular. 
It has been found convenient to select the following linear corrections and 
to arrange them in the order as listed: AVLp(Af 2# Af^, Af n , Af^ Af^, 


Af 1 


Af * 


A f 2 


AHifr ) ? and 


, Af^i); AVL H (AH*' n , AH^, AH X3 , ... 

sets of AVL^ ( AK£ n , AKfc 2 , AKfc^ , . AK kf 


AH T n ' 


AH, 


T W ' 


AfU 




2' 


AK]^ n ) . The nonlinear corrections are conveniently arranged as follows: 


AVNL_ 


Aa. 


Af 


Af 11 


Af', 

w 


Af l 


Af “ '); AVNL h (AH T v/ AH^ 1 , AHt 2 , AH» n _ i ) . and 


H # U4. n , ^ 

K sets of AVNL K (AK k ^, AK^, AK^, •• •# A^j^ . Here the nodal stations are 

sequenced from 1 at the wall (subscript w) to n at the outer edge of the 

boundary layer. The linear equations are conveniently sequenced as follows: 

L_ (linear boundary conditions and Taylor series expansions for f and its first 

and second derivatives); L u (linear boundary conditions and Taylor series ex- 

H. 

pansions for H T and H£) ? and K sets of 1^. (linear boundary conditions and Tay- 
lor series expansions for and I?£) . The nonlinear equations are sequenced as 
follows: NL„ (the momentum equation evaluated between each neighboring pair 
of nodal stations together with the two nonlinear boundary conditions and 
a„ constraint) ; NL U (the energy equation evaluated between each neighboring 
pair of nodal stations) ; and K sets of NL R (each of the elemental species 
equations evaluated between each neighboring pair of nodal stations) . The 
form of the resulting matrix equation is shown in Figure 3. Here, for exam- 
ple, [[ANLpJ and ^BNLpJ are matrices representing the coefficients of the 
corrections ^AVLj and ^AVNL^], respectively, arising from the nonlinear set 
of equations 
matrix ^ENLpj] . 

The first step in the matrix solution procedure is to invert the sub- 

and [AL ] 

cir r_r _ 

for p = F, H and K. 
p = F and H 


NL„ with the corresponding errors given by the single column 
F 


matrices £aL d J and to form the matrix products [AL ] -1 £bl ] L -— 


c el p h 


PP . 

p 


The former products have to be done only for 
since the linear equations relating the k th elemental species 
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AVNL. 
































to its derivatives (L R ) have the same form as the linear equations relating 
total enthalpy and its derivatives (L^) . Furthermore, this has to be done 
only once for a given problem as the matrices £aL ] and £bl J depend only 
upon the boundary layer Ti-spacing which can remain fixed as a consequence of 
the stretching parameter a 


H* 


The linear corrections can then be expressed in terms of the non- 
linear corrections [AVNL p ] and the linear errors as follows: 


r AVL 1 = - 

AL 1 

r bl 1 

I AVNL 

L p J 

PPJ 

L ppJ 

L pJ 

I 

Ixl 

lx J 

i 


r i - 1 

r q 

AL I 

-EL I 

L ppJ 

L p J 

Ixl 

I 


(151) 


where I = 3N - 2 and J = N + 3 for p = F and I = 2N and J = N for 
p * H or K with N the number of nodal points in the boundary layer. These 
can be introduced into the nonlinear equations to yield the reduced problem: 



r 1 


" 

BNL 

AVNL 

— 

ENL 


I xJ J I 


(152) 


where I = (K + 2) (N - 1) +3, 
matrix f BNL| are given by 


BNL ij 


= BNL ij " 


J = (K + 2) N + 3 , and the coefficients in the 

i [*vlf ltj _ ii53 > 


where p = F, m = 3N -2 and r = s = 0 for 1 ^ j ^ (N + 3); p = H, m = 2N, 

r = 3N -2 and s = N + 3 for (N + 4)^ j ^(2N + 3); and p = K, m = 2N, r = 3N - 2 + 

2NK and s = 3 + N (K + 1) for [(K + 1) N + 4]^ j ^ QK + 2) N + J] for the 

choice of linear and nonlinear corrections listed earlier. The coefficients 
in the £enlJ matrix are given by 

HE. . t *»i.l«|[«Vp] ' [-“pp]L 11541 

for p = F, H and K with the m and r having the same values as in Eq. 

(153) for each p. 
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Each time a coefficient in one of the original nonlinear equations (i.e., 
an ANL^ j or a BNL^) is formed, its contribution to [bNlJ and [enl] can 
be computed. Thus, the rather large matrix [ANlf] never has to be stored. In 
fact, it is highly significant that the only major blocks of coefficients 
which must be stored for the representation of all of the linear and nonlinear 
equations and for their solution are QbL fj J which is 3N - 2 by N + 3 , ^BL^] 
which is 2N by N, and [BNlf| which is [ (K + 2) (N - 1) + 3]] by £ (K + 2) N + 3] 
where N is the number of nodal points and K is one less than the number 
of elemental species. This should be contrasted with the size of the matrix 
of the complete set of linear and nonlinear equations which is Q(7 + 3K) N 4* 1 ] 
square . 

The matrix Eq. (142) is substantially reduced further as follows. First, 

the columns are rearranged so that the nonlinear corrections can be divided 

into two sets: AVNL a ( Aa R , Af “ , Af^', Af " Af^' ', AH*^, AHi^, 

AH" , Alt’’ , AKA' AKA 1 ) and AVNL, (Af , AH m and the aX v ) . Eq. (152) 

T n-1 k w k 2 n-1 b w T w k w 

can then be expressed as 



(155) 


where [_BNL a J i s a square matrix, being the coefficients of the I corrections 
[AVNL], with I = (K + 2) (N - 1) +3 and J = (K + 2) N + 3. Utilizing the 
same matrix reduction procedure employed previously (in going from Eq. (150) 
to Eq. (152)), the []AVNL a ] can be expressed in terms of the reduced set of 
corrections QavnLjJ as 


[ avm J - " 

BNLb] [AVNlJ 


I X x X 

I X J J 


+ [^a] 

enl! 

(156) 

Ixl 

. J 

i 


where I = (K + 2) (N - 1) +3 and i 

J = K + 2. 



The reduced set of nonlinear corrections [[AVNL^J ( Af w , AH Tw and the AKj^ 
are obtained from a consideration of the nonlinear wall boundary conditions. 


64 


Once these are determined, the remaining nonlinear corrections QavnlJ are 
obtained directly by use of Eq. (156) . The linear corrections ^AVL^J are 
then calculated using Eq. (151) . These linear and nonlinear corrections are 
then added to the corresponding primary variables in accordance with Eq. (115) , 

A_"U. 

thus completing the m iteration. The magnitude of the errors are checked 
and the procedure advances into the m+l^ iteration if the absolute errors 
exceed prescribed upper limits. If not, the iteration is completed for the 
current value of the streamwise position § and time t. 

Recurrence formulas for nonlinear wall boundary conditions are given by 
Eqs. (142) through (144) for assigned and by Eqs. (146) 

through (148) for steady-state carbon ablation, or can be obtained from any 
other consistent wall boundary condition. It is apparent that these equations 
must be expressed in terms of the [AVNI^]]. When P w v w (or f J) , H Tw (or T w ) , 
and the K kw are assigned, the ^AVNLjJ are given directly (e.g., Eqs. (142) 
through (144) ) . However, when fluxes are assigned or mass and/or energy bal- 
ances are required, the Ajft^ and/or Aq& w (given by Eqs. (128) and (126), 
respectively, evaluated at the wall) are functions of Qavl] and [AVNLj as 
well as [AVNLjJ. That is 


AFLUX i 


I 


3PLUX i 

5VLj 


AVL 


j 



9FLUX i 

avNL = 

3 j 


AVNL a 

a j 


+ 


c — > a FLUX. 

j 3 


(157) 


where FLUX. can be j* and/or q* and the summations are performed 
i K w a w 

over linear and nonlinear variables with nonzero coefficients. The [AVLJ 
and [AVNLj can be eliminated from these relations by making use of Eqs. (151) 
and (156). For this purpose it is convenient to look upon Eqs. (157) as 
additional nonlinear equations in the matrix Eq. (150) : 


[afluxI 

Faf ! 

bfJ 

AVL 

L J 

I 

L i 

X x J 

AVNL 




(157a) 


where I is the number of flux equations of concern, J = (7 + 3K) N + 1, 
AF ± j = BFLUX^/dVLj and BF^ = dFLUX ± /dVNLj . First, Eq. (151) is utilized to 
eliminate the [^AVL^ from Eq. (157) , yielding the result 


INFLUX 


I 



I X J 


AVNL 


J 



(158) 
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where I is the number of flux equations, J = (K + 2)N + 3 (the number of 

nonlinear variables), and the coefficients are given by 


BF ij 


hi 

■ “id- Z »i.i« [*vT K p]i 


(159) 


1= 1 


l, j-s 


EF i " " E t AF i.i + r Kp] _1 [- EL pp] 


(160) 


p Z-- 1 


where the subscripts are the same as those defined after Eqs. (153) and 
(154) . Next, Eq. (158) is rearranged as 


A FLUX = ^BF a ! BfJ 


X x J 


AVNL 

3 


- £efJ 


(161) 


Equation (156) is then used to eliminate the [AVNL ] from Eq. (161) , yielding 


[AFLUxJ = [BFj ^AVNlJ - J^EF^j 


(162) 


X x J 


where I is the number of flux equations, J = K + 2, and the coefficients 
are given by 

m 


BF b. . = BF b.. 

il 13 


, - & - J Kl " [-J L 


(163) 
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(164) 


with m = (K + 2) (N - 1) +3. For clarity, it can be seen that the matrix 
Eq, (162) is equivalent to 


66 


+ c k Af w + \ 


f 



- Ac *£ 


K \ 

E + V h t 

kk=l w w 
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E a K+l,kk A \k 

1 . 1 . W 


kk=l 


b K+l AH T + C K+l Af w + d K+l 
w 


> (162a) 


where the coefficients are constants during the m^* 1 iteration. Thus, cor- 
rections in the wall fluxes have been reduced to linear functions of correc- 
tions in wall state and total mass flux into the boundary layer (utilizing 

Eq. (145) and noting that Hm = h ) . 

x w w 

Equations (162) can be substituted into the recurrence formulas for the 
energy and mass balances (e.g. f Eqs. (146) and (147)) for steady-state carbon 
ablation. The resulting equations together with the additional recurrence 
formula (s) for the nonlinear wall boundary conditions yield the following 
matrix equation: 



I x I 




(165) 


where I = K + 2. The AVNL^ are then determined by matrix inversion 



C D 


1 x I I 


(166) 


It should be noted that the time required to invert this matrix is trivial. 
The time-consuming inversion is that of [bnlJ , required for Eqs. (156), 
(163) and (164), which has the dimension (K + 2) (N - 1) + 3. 

In this section the Newton-Raphson recurrence formulas for the linear 
and nonlinear equations have been developed. An efficient method for solving 
these equations has been outlined. The analysis is completed for the case of 
a boundary layer coupled to steady-state mass and energy balances. The pro- 
cedure utilized for coupling to a transient solution for internal conduction 
is described in Ref. 10. 


67 


SECTION 6 


RESULTS FOR INCOMPRESSIBLE AND COMPRESSIBLE 
SINGLE -COMPONENT BOUNDARY LAYERS 

The equations presented in Section 3 have been programmed in Fortran IV 
utilizing the numerical procedure described in Sections 4 and 5. Solutions 
have been compared to available results for several incompressible and com- 
pressible single-component boundary -layer problems as tests of the accuracy 
and convergence of the procedure. It has been seen that reasonably accurate 
results can be obtained with as few as five nodal points (the wall, three 
internal points, and the boundary-layer edge) and that three to four place 
accuracy can generally be obtained with seven points, although about 11 points 
have been required for some severe tests to obtain this level of accuracy. 
Convergence has been consistently satisfactory, four or five iterations being 
required for a starting solution and three iterations generally being adequate 
for subsequent (downstream) solutions. Some typical convergence and accuracy 
checks are presented in this section. 

Typical convergence of a velocity profile in an incompressible similar 
boundary layer with adverse pressure gradient is shown in Fig. 4 in terms of 
the conventional* Levy-Lees transverse coordinate. Starting with an assumed 
linear profile, the first iteration established the basic shape of the pro- 
file, the second iteration brought the solution within two to three percent, 
and the third iteration yielded results which were converged to three signi- 
ficant figures and compared favorably with the tabulated results of Loit- 

. . . 26 
sianskn. 

Validity checks for the nonsimilar incompressible boundary- layer problems 
of uniform blowing and uniform suction on a flat plate are presented in Fig. 5. 
The wall shear function, f", for these two cases is shown in Fig. 5(a) in terms 

" 2*7 

of the stream function at the wall, f .** The solution of Lew and Fanucci 

w 

is shown for comparison. Convergence required about three iterations at each 
streamwise station. Profiles of velocity ratio, f 1 , and shear function, f", 
are presented in Figs. 5(b) and 5(c). It is of interest that both blowing and 


*The normalizing parameter a H has been removed from r\ subsequent to the 
calculations in the accuracy checks presented in this section. 

**It should be noted that f v is zero at the leading edge of the plate, in- 
creasing with streamwise dimension for the case of suction and decreasing 
with streamwise dimension for the case of blowing. 


68 


VELOCITY RATIO, u/u 



Figure 4. Iteration History for a Similar Incompressible Boundary 
Layer Solution Near Separation (p = -0.19, Separation 
at -0.1988) 
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Profiles 


r 


suction results were efficiently obtained with the same r\ spacing, in spite 
of a ten- fold variation in boundary- layer thickness, by the use of the coor- 
dinate stretching parameter, a H , introduced in Section 3. 

The wall shear function for uniform blowing into an incompressible bound- 
ary layer is presented in Fig. 6 as a function of streamwise distance for a 
pressure of 1 atmosphere, temperature of 2000°R, and mass injection rate of 
0.005 lb/sec ft 8 . Results are also shown for the case where blowing is ter- 
minated at a streamwise position of 2 feet. It may be of interest to note 
that the nonsimilar effect of upstream transpiration decays quite rapidly, 
but that some influence persists for an appreciable distance downstream. 

Velocity profiles, shear function profiles, and maps of enthalpy ratio 
versus velocity ratio are presented in Fig. 7 for incompressible and compres- 
sible similar boundary layers with Prandtl number of unity and various posi- 
tive and negative pressure gradients and wall-to-edge enthalpy ratios. The 
results compare favorably with those of Hartree 28 and Cohen and Reshotko. 29 
These results were obtained with a 7-point nodal network. Two to five (gen- 
erally three) iterations were required for each problem, where the previous 
result was in each case employed as a first guess. 


Calculations were made for several single-component compressible boundary 
layers with variable properties corresponding to that of air at moderate tem- 


peratures 
of 




P °.85 


and 


C « T 

u e V2H T 


o.l 9 


\i « T 0 * 70 and Pr^ = 0.7) for various values 
As an example, profiles of velocity ratio, 

for P = 1 atm, 

0 for several values of u s /2H m . Conver- 
- - e T 

gence was similar tS that obtained in the Cohen and Roshotko compfrisons. 


and Pr = 0.7) 

H T w / ' H T e 

shear function, and temperature are presented in Fig. 8 


T w = 1200° R, H t /H t = 0.2 and 3 = 


As an example of a nonsimilar compressible boundary-layer solution, re- 
sults are compared in Fig. 9 to results obtained at the NASA Ames Research 
Laboratory with the Smith and Clutter finite-difference procedure 8 for the 
problem of Mach 10.4 flow over a 7^ degree cone with uniform injection down- 
stream of the 0.1574 foot station. The wall shear is compared in Fig. 9(a) 
whereas representative velocity profiles are compared in Fig. 9(b). The re- 
sults agree to nearly four significant figures for the similar no-blowing 
solution and the first few stations downstream of the point where mass injec- 
tion is initiated. Further downstream, the results are still in reasonably 
good agreement, considering the ^-derivatives and possibly the f^ and ? 
are computed somewhat differently by the two methods. The present results 
were invariant with streamwise spacing. The effect of nodal distribution 
across the boundary layer was not investigated. The effect of streamwise step 
size on Smith and Clutter results was not available. The first solution re- 
quired 4 iterations, whereas the downstream stations, including those near the 
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Figure 6. Effect on Wall Shear of Upstream Air-to-Air Injection on a Flat Plate 





TRANSVERSE COORDINATE 



(a) Velocity Profiles 

Figure 7. Solutions for Similar Compressible Boundary Layers 
with Prandtl Number of Unity for Several Values of 
P and Wall-to-Edge Enthalpy Ratios 
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SHEAR FUNCTION. 


V 



(b) Shear Function Profiles 
Figure 7. Continued 
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Temperature 



Figure 8. Profiles 
Temperal 


Layers 
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Figure 9. Air-to-Air Transpiration into a Compressible Boundary Layer over a 7.5° Cone 
with an Impermeable Sharp Nose Tip 









zero wall-shear condition, required three iterations each. The entire solution 
was generated in 45 seconds on the IBM 7094 computer, an order of magnitude 
less than the results obtained with the Smith and Clutter procedure. The 
stability of the present procedure is indicated by the fact that no difficulty 
was incurred in obtaining solutions down to the zero-shear condition. 

In order to gain some confidence in the ability of the solution procedure 
to compute physical boundary-layer thickness, results were correlated with 
experimental data reported by Watson et al^ for laminar shock interaction 
over a flat surface. A velocity profile is presented in Fig. 10. The correla- 
tion is good with the exception of values near the wall where it is reasonable 
to expect some experimental error due to probe-shock interaction. 

Results were obtained for a boundary layer in a ficticious gas with rap- 
idly varying properties (e.g., Prandtl number variation from 0.6 to 10 across 
the boundary layer) in order to study convergence in an extreme situation. 

A converged solution was obtained in 7 iterations, starting with the usual, 
uninspired, built-in first guesses. 

The problem shown in Fig. 8 (where Mach number effects were investigated) 
with moderately varying properties and the problem of rapidly varying proper- 
ties were repeated with various partial derivatives which are used in the 
iteration process but do not appear in the boundary-layer equations (such as 
dPr/dh and dC/dh) set equal to zero. In the former case, the total number 
of iterations was increased by two for the entire set of five solutions. In 
the more severe case, 11 iterations were required (whereas 7 were required 
previously) . The final results, of course, were unaffected. 

SECTION 7 

SOME RESULTS FOR MULTICOMPONENT BOUNDARY LAYERS 

To date no accuracy studies have been performed for multicomponent 
chemically-reacting boundary- layer problems. However, convergence has con- 
sistently been satisfactory. To illustrate, an iteration history is given in 
Table IV for graphite ablation in air. In this problem, the ablation rate was 
assigned and surface temperature was determined by a coupled mass balance 
at the surface together with heterogeneous equilibrium. Considering that 
the initial guesses for velocity ratio, total enthalpy, and elemental mass 
fractions are built-in, linear profiles with respect to r), and that there 
are no constraints on the size of the corrections, convergence such as that 
shown in Table IV is highly encouraging. Furthermore, it is significant that 
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DISTANCE NORMAL TO PLATE (IN) 



Figure 10. Laminar Interaction on a Flat Surface: = .10.5 T 



this convergence has been obtained while ignoring all second derivatives and 
some first derivatives in the iteration process. 


TABLE IV 

CONVERGENCE FOR A TYPICAL GRAPHITE -IN- AIR 
ABLATION PROBLEM 

(Assigned Ablation Rate, Surface Conditions 
Determined by Coupled Mass Balance) 


Maximum Relative Errors* in 


Iteration 

Normalized 

Wall 

Shear 

Wall 

Enthalpy, 

Btu/Lb 

Wall 

K " p r”' 

Momentum 

Equation 

Energy 

Equation 

Elemental 

Species 

Equations 

1 

0.3135 

-223.0 

1992 

1.0 

3.7 

0.31 

2 

0.1727 

-804.0 

1469 

0.29 

0.77 

0.13 

3 

0.1628 

-795.2 

1478 

0.10 

0.53 

0.025 

4 

0.1670 

-794.2 

1479 

0.0064 

0.0028 

2.7 x lO"" 

5 

0.1670 

-794.3 

1479 

6.9 x lO -4 

1.9 x 10“ 3 

1.3 x 10" 1 


* A relative error of 1 x lO”* 3 corresponds to nominal 4-place accuracy 


Boundary- layer profiles of velocity ratio, temperature, and elemental 
mass fractions are presented in Table V for an air boundary layer over a flat 
plate with unequal diffusion coefficients for all species, with and without 
thermal diffusion. Mole fractions are presented graphically in Fig. 11. For 
assumed equal diffusion coefficients and in the absence of thermal diffusion, 
the elemental mass fractions remain constant across the boundary layer at the 
assigned edge values of 0.23 and 0.77 for oxygen and nitrogen, respectively. 
Consideration of unequal diffusion coefficients for all species is thus seen 
to have a substantial effect on elemental mass fractions. With unequal dif- 
fusion, the elemental mass fraction of oxygen first decreases slightly and 
then rises to a maximum value of 0.2723 at the wall. When thermal diffusion 
is also taken into consideration, the elemental mass fraction at the wall is 
decreased slightly to a wall value of 0.2709. These wall values of elemental 
mass fractions necessary to maintain zero mass flux at the wall are analogous 
to the adiabatic wall temperature for zero heat flux at the wall. 
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TABLE V 


AIR BOUNDARY LAYER OVER A FLAT PLATE 
WITH UNEQUAL DIFFUSION COEFFICIENTS: 


(p = 

= 1 atm, H t 

= 5000 Btu/lb, 

T = 5000° 
w 

R, s = 1 

foot) 

r\ 

y. ft 

u/ u e 

K o 

■*N 

T,°R 

Thermal 

diffusion neglected 




0 

0 

0 

0.2723 

0.7277 

4000 

0.584 

.00075 

0.1889 

0.2626 

0.7374 

5355 

1.168 

.00171 

0.4044 

0.2475 

0.7525 

6165 

1.752 

.00284 

0.6148 

0.2308 

0.7692 

6861 

2.336 

.00418 

0.8000 

0.2239 

0.7761 

7722 

3.504 

.00735 

0.9699 

0.2269 

0.7731 

8892 

5.840 

.01427 

1.0000 

0.2300 

0.7700 

9204 

Thermal 

diffusion included 




0 

0 

0 

0.2709 

0.7291 

4000 

0.584 

.00075 

0.1887 

0.2594 

0.7406 

5350 

1.168 

.00171 

0.4040 

0.2447 

0.7553 

6170 

1.752 

.00284 

0.6145 

0.2288 

0.7712 

6877 

2.336 

.00418 

0.8000 

0.2240 

0.7760 

7748 

3.504 

.00736 

0.9700 

0.2279 

0.7721 

8903 

5.839 

.01429 

1.0000 

0.2300 

0.7700 

9204 


This effect could be significant in ablation problems since substantially 
more oxygen is available at the wall for reaction. Therefore, it is pertinent 
to investigate the cause of the observed behavior. As a consequence of the 
lower dissociation temperature of 0 2 relative to N 2 , the oxygen is almost 
completely dissociated into the relatively more mobile atomic oxygen at the 
edge of the boundary layer, whereas the nitrogen is only slightly dissociated 
(see Fig. 11) . Moving from the boundary- layer edge toward the wall, the ni- 
trogen recombines more readily than the oxygen. As a consequence, the gradi- 
ent in the atomic nitrogen concentration, although small, is greater than 
that for atomic oxygen until an r\ of 2.5 or so is reached. Thus there is 
a small net flux of nitrogen inward from the boundary-layer edge with a con- 
sequential small decrease in oxygen elemental composition. Closer to the 
wall, substantial recombination of the oxygen occurs, producing a large gradi- 
ent in atomic oxygen concentration and thus a flux of oxygen toward the wall. 
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Thermal diffusion produces a small effect in the present problem since 
temperature gradients are relatively small. Calculations performed at high 
edge enthalpies have shown large thermal diffusion effects, as a consequence 
of the larger temperature gradients, resulting in wall oxygen concentrations 
below edge values. 

Profiles of velocity ratio, temperature, shear function, and mole frac- 
tions across a boundary layer into which a large quantity of Apollo heat- 
shield material is being injected are presented in Fig. 12. These results 
were obtained for an assigned surface temperature and assigned component 
fluxes (irig and m c ) and utilized a 30-component chemical model. A converged 
solution was obtained in 7 iterations, starting with an air boundary-layer 
solution with the same wall temperature and same edge conditions but with 
no mass injection. The convergence histories of the wall shear function and 
maximum relative errors are presented in Table VI. In this calculation, the 
corrections in the elemental species equations were not allowed to exceed 0.30. 

TABLE VI 

CONVERGENCE FOR ABLATION OF THE APOLLO HEAT-SHIELD 
MATERIAL INTO AIR 

(Assigned Surface Temperature and Component Fluxes) 




Maximum 

Relative Errors* in 

Iteration 

Normalized 

Wall 

Shear 

Momentum 

Equation 

Energy 

Equation 

Elemental 

Species 

Equations 

First Guess 

.2007 

— 

— 

— 

1 

.0883 

0.26 

0.16 

0.30 

2 

.1343 

0.19 

0.42 

0.30 

3 

.1016 

0.095 

0.38 

0.23 

4 

.0995 

0.069 

0.21 

0.11 

5 

.0987 

1.4 x 10" 3 

C .045 

0.026 

6 

.0987 

4.8 x 10“ 6 

6.5 x 10“ 3 

5.1 x 10 -3 

7 

.0987 

2.6 x 10 -6 

9.0 x 10"* 

7.5 x 10 -4 

*A relative 
4-place acc 

error of 1 x 10~ 3 corresponds to nominal 
curacy 


The nonsimilar boundary layer around a sphere-cone reentry body with 
water injection was studied to determine the injection rates required around 
the body to maintain a uniform wall temperature of 5000° R. A 16-component 
chemical model was employed in these calculations. The distribution of water 
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Distance from surface, inches 
(a) Velocity Ratio, Temperature, and Shear Function 

Figure 12. Boundary-Layer Profiles Over the Apollo Heat- 

Shield Material: Assigned Wall Temperature and 

Component Mass Fluxes (m^ and m c ) 
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Mole fractions 









s 


injection rates and the distribution of velocity profiles around the body are 
presented in Figs. 13(a) and (b) , respectively. These solutions, which in- 
cluded evaluation of edge conditions, a similar solution at the stagnation 
point, and nonsimilar solutions at ten additional stations, were obtained in 
a PP rox i ma tely 2.5 minutes on an IBM 7094 computer. Approximately 60 percent 
of this time was spent in the equilibrium chemistry subroutines, it is thus 
pertinent to mention that the total computational time could be substantially 
reduced by the use of a specialized chemical procedure for the particular sys- 
tem of interest. Also, substantial time was spent in tape operation intro- 
duced by the use of overlays . 
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Figure 13 
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APPENDIX A 


INTRODUCTION OF THE APPROXIMATION FOR MULTICOMPONENT 
THERMAL DIFFUSION COEFFICIENTS INTO 
DIFFUSIVE FLUX RELATIONS 


In this appendix the approximation for multicomponent thermal diffusion 
coefficients is introduced into the expression for diffusive mass flux of 
species i, diffusive mass flux of element k, and diffusive heat flux. 

The diffusive flux of species i in a multicomponent gas incorporating 
the bifurcation approximation for binary diffusion coefficients can be ex- 
pressed as Eq. (23) : 


PDH 2 

V7W 


az. 

9T + {Z i ' K i } 


a in n 2 

ay 


+ 


a in T 
D i *y 


(A-l) 


when the diffusion factors, F i# are considered to be invariant with tempera- 
ture. Introducing the approximation for multicomponent thermal diffusion CO' 
efficients (Eq. (26)) 



c t pD U 2 


< Z i 


K i> 


into Eq. (A-l) yields 


p^ 2 

-azi 

K - k i) 

a in n 2 a in T\" 

3 i “ 

dy 

ay 1 c t ay / 


(A- 2) 


(A-3) 


or, equivalently. 


pdp- faz. 

3y + ^ Z i 


K i) 


c t 

a in (m 2 t ) 

ay 


(A- 4) 


which is the desired form for the diffusive flux of species i including the 
approximation for unequal diffusion coefficients embodied in Eq. (26) . 
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Likewise, the diffusive mass flux of element k is given by 


where 


k Wtfi |^3y + (2 k V 3y + D k 




Substituting Eq. (A- 2) into Eq. (A-6) produces 


S In T 

ay 


C.pDn, 




a ki (Z i - K i> 


or by the definitions of Eqs. (16) and (25a) 


T C t pI *2 


i 7JT (z k - V 


The elemental diffusive flux- can thus be expressed as 


i * , , L. 

■’k Ujflj ay + (Z k V 3y 


a in (n 2 ^ ) 
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The diffusive heat flux can be written as 


< 3 =. = 


<«. * V, 


"M T v; ay 


+ <* + pe H C p ) |2- + 




h i - S > (a- 10) 


where 


S s 




^D, j. j. 


k d 


i D 


Substituting for D ± and the term S becomes 


S = 




b .it 

K i K i 


i 3 


(A- 11) 


(A-12) 


Replacing the mole fraction x_. by and utilizing the definition 


Z i 3 ^K/ F iti 2 » s becomes 


c t RTp 2 r- F L , ^ 

= *1 L L * j \ |^ F i 


i 3 


'j. - — j • 

J i Kj ^3 


(A-13) 


Expanding 


S = 


c t RTM 2 


I VV !hfi ii _ V V 
^ 2 Z. ^ \ L L n j % 


1 3 


1 3 


»i ' LL »i ») J 


1 3 


i 3 


(A- 14) 
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* The second and fourth terms are identical upon interchange of subscripts but 
of opposite sign and hence cancel each other. Utilizing 1 / 7/1 = ^ K _■/???_■ 
and the definition of = ]^ x i F i' S bec onies 


S 


V i y j i F i 

Ml V ^ m \ 


Rearranging 


S 



(A-15) 


(A- 16) 


Substituting Eqs. (A-4) and (A-16) into Eq. (A-10) , the expression for 
diffusive heat flux becomes 


% - i p(e M + v) t£ uV2) + + If 


r -1 9K- pDp v-> faz. 

+ pe D h i ay - + Z. + (Z i ' K±) 


d jen (h 2 T t ) 

3y 


Now 


X 



c RT 

+ U^- 


(4l 




Z 


a (K.h.) 
ay 


z 


K. 


1 




= . r IS 

ay p 3y 


(A- 17) 


(A- 18) 
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Utilizing this and eliminating all summations in favor of defined system 
quantities yields the following result for q & 


q a = * 


|> (e M + 


v) 


3(u 2 /2) 

Ty 


+ P e ^ It7 - C |S + 


\ p*i 2 

il _ I £ + c t 3R | 

1 ^ 

Jy \ p W 2 I 


— - ^ T 

+ + ay 


ftT 9,j 3 

If + 't™ W 


d in (n_T r ) 
+ (h - h + c t RTn 3 ) ^ 


(A 


where \x 2 = £ ^i x i/ F i» U 3 = £ h = and C p - E Z i C p i* 




-19) 
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APPENDIX B 


DERIVATION OF THE TRANSFORMED NONSIMILAR 
LAMINAR BOUNDARY-LAYER EQUATIONS 


In this appendix, transformation relationships for the Levy-Lees trans- 
formation, modified by the use of a stretching parameter, are developed and 
applied to the expressions for diffusive fluxes and to the boundary-layer 
conservation equations. 


The desired independent variables J and j) are defined by Eq. (33) 


T = e 


where a„ is a function of £ only, and 


n = a- 

a H 


& 

S =J u e Pe^e r o K ds 


k y 
r o u e r 

(2S) h Jn 


p dy 


(B-l) 


(B-2) 


( B-3) 


In addition, f is defined as Eq. (35) 


-f s f a. 
W L U e 


d n 


(B-4) 


from which it follows that f' = u/u g where the prime refers to partial 
differentiation with respect to r\ . 

The old partial derivatives can be expressed in terms of the new partials 


by 


a 

_ a 
y' 3 ? 

¥ 

n * 

y + ^ 

t 

a 

_ a 

ae 

, a 

dn I 

3y 


n * 

s + ^ 

e 


where 


If 


If y = PeWl* 


If 


K 

Vo' 

y l» ' ue? 5 " 


(B-5) 

(B-6) 

(B-7) 

( B-8) 
(B-9) 
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Furthermore 


where 


such that 



a _ s_ | f^H j 

n df - a H d£ I dr) 1 

S _ d_ 

^ 5 "-2b 3n|-j— 



( B-10) 

( B— 11) 

(B-12) 

(B-13) 
(B-14) 
( B-15) 

(B-16) 

(B-17) 

(B-18) 
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The global mass conservation equation (Eq 


( 3 ) ) can be integrated to yield 


pv 


1 d / H C . 

- Pw v w = - 7 SS r o I P u dy 

r o 


Differentiating r) for constant § , introducing f ' , solving for 
(2£) 2 f* dq , and integrating yields the result: 


( 2 ?) 


' 2 j f dt) = r oJ~ p u d y 


Utilizing this and the definition of f , Equation (B-19) can be written 


P V P W V W " r K "§S 

o 


( 24 ) 15 ( f - f w ) 


With Eqs. (B-5) and (B-8), Eq. (21) becomes 


pv - p v v w - - v] - Is | 


Performing the differentiation 


P v P w v w p e U e^e r o 


(2 4) ^ 

L (2?) 2 


d£ df w 

Tt ~ W 


< 2 ^ li 


For convenience, define 


PeWo 

(26)^ 


df 

f + 2£ ~ 
w s d£ 


Then 


Ji r 


Pe u e^e r o 

pv r — 

(2?) 2 


£ ♦ (26) |f 


T & 


B-19) 


B-20) 


as 

B-21) 


(B-22) 


( B-23) 


(B— 24) 


(B-25) 
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Utilizing this, plus Equations (B-5) through (B-9) 


5 , d _ 2 2%- c i ^ [ f . 5 

pu 31 pv ^ " °Pe u e^e r o f 3? “ I 2f + 3? 3^ 

y s L rj ' 


( B-26) 


which is the desired form for the operator. 

In accordance with Equation (B-24) , the use of this operator requires 

that f be defined as 
w 


ju r p w v w 

f w = - < 2 *> / " X 

J o Pe u e^e r o 


(B-27) 


The diffusive flux of species i is given by Eq. (27) 


pD p- [to, b In Ip., T 

>i = - H^TT 3F + (7 'i - K i> T7 


(B-28) 


Applying the transformation of Eq. (B-6) with the aid of Eqs. (B-7) and (B-9) 
yields 


Pe U e^e r o C C t 


( B-29) 


where Sc is defined as \l^\ 3 lW pD|r 2 and the prime refers to differentiation 
with respect to rj . Utilizing Eq. (B-17) : 


afp W * '** ' Kl, ['”(“ 2 T#l l 


(B-30) 


where the prime now refers to differentiation v/ith respect to r\ . 

The laminar form of the streamwise momentum equation is given by Eq. ( 5 ) 


du , du _ d Sul dP 

P u SS + pv 3? * 3y 3?J “ 'Si 


(B-31) 
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Applying Eqs. (B-26) , (B-5) , and (B-6) yields by direct substitution: 


2 2K 

* , du 

f 4. 1 

du 

u e r o P 5 

u e r o P du 

U e ^e r o 

f s? " 

2? + 3?) 





„ , 2 k Sp Sn 
“ P e u e^e r o *5% " 


Sp 


(B-32) 


Now u “ f ’u where u is a function of £ only 
e e ^ J 


= u |£1 = u f 
ot| e drj e 


( B— 33) 


du _ ,, du e , „ b£< 

3?" f W U e W 


(B-34) 


So that 


2 2K 

PP e u | P e r o 


/•(*•£♦ v] 


u e r o P 


V u e r o P u e f " 

3 


(2|) 


“ Pe u e tI e r o H If " !• y P ' 


(B-35) 


where the prime denotes differentiation with respect to r\ . Multiplying by 

( 2 ?)/pP e u e M e r o K y ields : 


( 2 y f,a Jr + ( 2 £) fl If 1 “ f f " “ ( 2 ?) f " |f 


-OU-- f M - 121L 4^ - 

DU 2 dj 

Pe^e / pu 


(2?) 


3 2 * 

PPe u e ^e r o 


P' (B-36) 


Note that all properties evaluated at "e" are functions of £ only. 
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and 3 = 2 


yields 


Defining C = pn/p \l 


i_ du e 

^ W 


(2£) 

(Cf")' + ff" - 9f ,a - 1211 4 t 

2 O c 3 oft 

P u e PPe u e ^e r o 


^7 P 1 = (2?) 


f 41^- - 41 f 
- W 3f f J 


The final term on the left-hand side of this equation involves 4-; 
can be expressed with the aid of Eq. (B-8) as 


ZK 

= p u ll r 

*e e^e o 


* 1 , 


Therefore 


(Cf") ’ + ff" - 3f ,a - liii 

P u « 




(B-37) 

which 

( B-38) 


■ lr - If £- ] 

Introducing the additional transformation of Eqs. (B-16) and (B-17) 


1_ 

o K. 


C ^ f" 

■ 

+ f 


- 3 

\ f,2 l 

a 

L H J 


Qt 

1 H ' 


CL 

H 1 


— da T 


.p.i.ZL + p, a 

p 1 ■n’T' r>__ — 


p u . 


d? “h d§ 


d_ / r o u e r , 

J. P y 


~ f ‘ 

' — (~ f •! 

_ / 

( - f| 

1 — 

] a H 

L dfK 1 

1 a H df , 

l l a H I 



df _ 

' a_ 

da H 

1 -F ■ 

1_ f" 

bl ~ 

1 a H 

d? 

r 

— ! 

7 

a H 

> 


(B-39) 


(B-40) 


where the prime denotes differentiation with respect to t] and use has been 
made of the relation 
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(B-41) 


Hr = °» I 


£ 

= a.. 

■s / 


H 

y 

W 


-% fpa y ]l 
° h < 2 ^ )J 


Multiplying by ex* and rearranging terms yields (noting that a R = a (£) 
only) : 


rcfii 

' 

) Sp _ 

JL. + P . 

's J 


% f . dy |1 

L a H J 

P< 

1* 

“h dl 

w 

a H 

/J 




= 2f 


f. Ml fu Sf f' 5 da H 

SI si a H dl J 


(B-42) 


which is Eq. (52) of this report. 

The normal momentum equation is given by Eq. (7) 


M - = o 

Sy r c 


(B-43) 


Applying Eq. (B-6) with the aid of Eqs. (B-7) and (B-9) yields directly 


u r p 2 

■ ■ ■ P’ - -fiH- = o ( B-44) 

(2l)^ r c 

where the prime denotes differentiation with respect to r\ Utilizing 
Eq. (17) and f' - a H u / u e 


P 


u f ■ 2 ( 2l) ** 


a H r c r o 


= 0 


(B-45) 


where the prime refers to differentiation with respect to t) . This is 
Eq. (54) of this report. 

In the event that normal pressure gradient can be neglected, Eq. (B-45) 
is replaced by p' = 0 and the compressible Bernoulli equation yields 



d? 


Pe 



( B-46) 
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Applying this to Eq. (B-42) , the streamwise momentum equation becomes 


ff" + 


+ 3 a* •— - f 1 2 


2 j\f< ill _ f if _ ill 

Sf 3? a H df 


which is Eq. (55) of this report. 

The laminar form of the energy equation is given by Eq. (30) 

ou + oV *r = 3 | a S(u a /2) + , St 


P u -gr + P v - -Sy 


. ~ . c t R &r . DIP ^3 

+ -Sr - \ c P + ~ 2 ) 'Sy + c t RT W 


+ (fi - h + c t RTn 3 ) 


a in l|J.-T 


Applying Eqs. (B-6) and (B-26) 


PP e u e^e r T f ’ -5T " (if + If) “t 


K K 

u e r oP u e r oP 


(2£P (2?) 


LI U 2 f‘f" 

^ e 


+ A T' 


p D U 2 I~ c a R 

M-l™ l 1 P HjPj 


T' + c t RT 


r i Min 

+ (h - h + c t RT p. 3 ) in p 2 T ^ I + q r 


where the prime denotes differentiation with respect to rj Multiplying 

by 2 ?/ op e u| M- e r o* and rearranging 
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( 




C ir ff" + 
e 


C C P T> , C 


Pr 


Sc 


c*R 

U- - c + 

P 


T ' + c^RT[JU 


\ I 


(h - h + c t RTp- 3 ) Jjenj p 2 T t 


q r (2gP 

Pe u e^e r o 


= 2? 


£ . 2* Hl^l 


where Pr = C p p/A and Sc = p.^ W pD p 2 . Applying Eqs. (B-16) and (I 
and multiplying by a H yields 


c 

f ' f" 

a H 

2 

L a H 


C 

a 2 + _E 
e Pr 


— — fh' - (c + - - |t' + c.RTpl 

aS \ 1 p ^1^2 j 3 


+ (h 


- h . o t RT. 3 ) [x»(./‘j] JJ + 


p U Li _ 

r e e^e o 


21 


Sh t 
f' - 

L H 


l a H dl / ^ ’ S? 


f- daj 

- 


f 1 

\H dT J 



2? 


f ~ - 14 . -rr 

L a? ^ d«. 


where the prime denotes differentiation with respect to t) . This is 
Eq. (56) of this report. 

The laminar form of the elemental species equation is given by Eq. 


^ 5l <k + v d «k _ d J ^2 
p -xt p w w - 


'SZt „ _ a in ( p 2 T 


l+ «k 


(B-50) 

-17) 


(B-51) 

(29) 

(B-52) 
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(B-55) 


where the prime denotes differentiation with respect to T) . This is Eq. (58) 
of this report. 
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APPENDIX C 


SPECIAL STAGNATION POINT CONSIDERATIONS 


In this appendix equations are developed which are applicable at the 

stagnation point of blunt bodies for the calculation of the flux normalizing 

parameter, a*, the streamwise pressure-gradient parameter, ft, and the wall 

stream function, f . 

w 

The flux normalizing parameter, a*, is defined by Eq. (44) as 


a* 


p p u r 
e'e e o 


( 2 §) 


(C-l) 


In order to compute the value of a* at the stagnation point, consider the 
streamwise pressure-gradient parameter, p, defined by Eq. (53) 


rr dU 

3 = 2 

u e d§ 


(C-2) 


The transformed streamwise parameter, §, is defined by Eq. (31) 


s 

? = / Pe u 


e^e r o 2Kds 


(C-3) 


Differentiating Eq. (C-3) and introducing the result into Eq. (C-2) yields 

du 


P = 2 


1 e 


u e ds p u p r 3H 
e^e o 


(C-4) 


Utilizing Eq. (C-l) 


P = 


du 


a* 2 


ds 


(C-5) 


Solving for a* yields the result 


a* = 


du e / 
p e^e ds7 P 


(C — 6) 
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The du e /ds| o needed for evaluation of a* if Eq. (C-6) is to be em- 
ployed can be obtained as follows if Newtonian flow is assumed in the vicin- 
ity of the stagnation point. In the case of Newtonian flow, 

p = P cos 3 6 (C-7) 

e e o 


Near the stagnation point sin 0 = s/R gff where R 0ff is an effective nose 
radius of the body. Hence 


e o 


'i.-i 


\ “eff 


(C-8) 


Utilizing the Bernoulli equation 


P o = P e + ?e T 


(C-9) 


Equation (C-8) becomes 


P a u 


e. e 


= P 


° R 3 
eff 


(C-10) 


or 


2P 


eff p e 


(C-ll) 


The density is nearly constant in the vicinity of the stagnation point. 
Hence 


du e 

1 

[ 2P °I 

ds 

~ -D 

o R eff 

l p 0 ) 


(C-12) 


In order to compute the stagnation-point 8, consider the definition 
of § (Eq. (C-3) ) in the form 
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(C-13) 



For Newtonian flow, comparison of Eqs. (C-ll) and (C-12) show that u e / s -► 
du e /ds in the vicinity of the stagnation point. In addition, r Q /s -► 1 
and pa is a constant. Hence 

G G 


du 


2k + 2 p e M e ds 


e s 3 *+s 


(C-14) 


Substituting into Eq. (C-4) yields 


P 


o 


1 

K + 1 


jl£M S 

u e ds 1 



1 

H + 1 


(C— 15) 


Hence, is unity for planar blunt bodies and 1/2 for axisymmetric blunt 

bodies for assumed Newtonian flow. 

The wall stream function, f w , is defined by Eq. (60) as 


w 



p v df 
K w w Z 


p u \JL r 
r e e e o 


(C-16) 


Differentiating Eq. (C-3) and introducing the result into Eq. (C-16) yields 


f„ - - (25) 




/ p v r K ds 
M w w o 


(C-17) 


This can be written as 


f = “ (2§) 


-H 


w 


s 

/ P V 
w w 


s*ds 


= - ( 2 §) 


-h 


s 

I 


_H+1 


— ; — t p v 
K + 1 w w 


ds 


K+l 


(C-18) 
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where again r Q /s ^ 1 in the stagnation region. The integrand in Eq. (C-18) 
is well-behaved, starting with a finite value at s = 0. However, at the 
stagnation point, the value of the integral is zero and f^ is indeterminant 
since § = 0. Thus, it is necessary to develop a special relation for f 
Applying l'Hospital's rule to Eq. (C-17) , ° 


w 


p v r K ( 2§) ^ 
H w w o v 

as. 

ds 


(C-19) 


Differentiating Eq. (C-3) and introducing the result into Eq. (C-19) 



o 


( 25 )^ 

Pe u e^e r o* 


P 


V 

w w 


Introducing Eq. (C-l) , 


= - 


(C-20) 


In summary, Eqs . (C-6) , (C-12) , (C-13) and (C-20) can be used to calculate, 

respectively, the a*, du^/ds, p and f^ at the stagnation point of a 
blunt body, and Eqs. (c-13) and (C-18) are useful for evaluating the integrals 
§ and f in the vicinity of the nose region. 
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APPENDIX D 


ALTERATION OF THE BOUNDARY-LAYER EQUATIONS TO ENABLE 
CONSIDERATION OF GENERALIZED BOUNDARY- LAYER- EDGE CONDITIONS 


In the conventional treatment of the hypersonic boundary layer, the 
boundary- layer-edge conditions are obtained from an inviscid solution as 
those conditions which exist on the "stagnation streamline" (i.e., that 
streamline which crosses the shock wave such as to become the stagnation 
point and wall streamline) . This same procedure can still be utilized to 
include boundary- layer-displacement effects if a new body shape is consid- 
ered according to the effective displacement of the flow due to the presence 
of the boundary layer. However, when entropy layer or nonadiabatic flow- 
field effects occur, this conventional approach is no longer adequate since 
the edge boundary conditions become functions of the local stream function 
as well as the streamwise coordinate. Hence, in these cases it is necessary 
to normalize the dependent parameters of the boundary layer in such a fashion 
that asymptotic solutions are achieved at the boundary- layer edge. A pro- 
cedure for accomplishing this task is presented in this appendix. 

The procedure used herein consists of normalizing the boundary-layer 
equations with respect to the f = 0 streamline when performing the Levy- 
Lees transformation described in Section 3 of this report. The equations 
then remain unchanged except that in all equations the subscript "e" 
is considered as this reference conditions, and "edge" as the actual edge 
condition. For example, the definition of f' becomes 


f ' = a HlT 

e 


= a T 


conditions are given by 


edge 


u I 

u 

u edge 1 

u 

e 

^edge 

is 

^edge* 

Cc 

p edge 

= 


(D-l) 


edge 


edge 


JL 

P 0 


(D-2) 


edge 


and so on, where u g , , M e# etc., are supplied by the inviscid solution 

for the f = 0 streamline, and ( u / u e ) e( ^g e and (p/P e ^ ec jg e are a ^ so supplied 
by the inviscid solution, but as functions of § and f. That is 
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In order to illustrate this procedure, consider the following sketches 
which represent typical profiles of ( u A e ) edge , u/u edge and f 1 at a given 
5 for two selected values of r l ed g e in the presence of an entropy layer: 



3 4 3 4 3 4 


^edge ^ ^ 

Sketch (a) Sketch (b) Sketch (c) 


In this example, the viscous effects are considered to be confined to an r\ 
less than 3 so that both h ed g e considered are out of the boundary layer but 
in the entropy layer. In Sketch (a), the ( u / u e ) edge is shown to increase 
with distance from the surface, as would be expected. Since u edge has dif- 
ferent values depending on the choice of rf . u/u . profiles would dif- 

edge edge r 

fer accordingly, as shown in Sketch (b) . The u/u e is# of course, indepen- 
dent of T l ed g e since u q is a given value for a given 5; however, as 
shown in Sketch (c) , the f'/& H does not approach unity at the edge of the 
boundary layer, and, in fact, depends on the value of h ed g e chosen to rep- 
resent the edge of the boundary layer. 


In the absence of an entropy layer or other such effect, u edge = u e and 
f’/a H does approach unity at the edge of the boundary layer as in the con- 
ventional solutions. It is for this reason that the f = 0 streamline is a 
convenient choice for the reference condition. 
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APPENDIX E 


ONE-DIMENSIONAL RADIANT HEAT FLUX IN AN ABSORBING 
BOUNDARY LAYER WITH ANGULAR- DEPENDENT INCIDENT RADIATION 


In this Appendix an expression is derived for the net one-dimensional 
radiant heat flux in the boundary layer normal to the surface, q r , which is 
needed in the energy equation (Eqs. (8), (30) or (56)). The wall value, 

q rw , which appears in the surface energy balance (Eq. (62)) is also presented. 
The multicomponent gas in the boundary layer is allowed to absorb diffusely 
as well as to emit, but scattering is neglected. The wall is assumed to 
emit spectrally and to reflect diffusely, but transmission of radiant energy 
is neglected. The incident radiation at the boundary layer edge is allowed 
to have an angular dependence to approximate intense radiation from the stag- 
nation region of the inviscid flow field or from a nuclear explosion outside 
of the boundary layer. The derived relations are an extension of the basic 
relations of Goulard and Goulard (Eqs. (17) and (18) of Ref. E-l) to include 
the specific boundary conditions described above. The present result also 
reduces to that of Cess (Eq. (5) of Ref. E-2) if the radiation layer is con- 
sidered to extend to infinity with no incident flux at this edge boundary 
condition. 

The one-dimensional approximation implies that the flow field extends 
uniformly to infinity on planes which parallel the plane tangent to the sur- 
face at the streamwise position of interest, but that properties may vary 
from plane to plane. It also means that the net radiation transfer along a 
ray which lies in one of these planes and passes over the streamwise posi- 
tion of interest is constant but may differ from the value along a different 
ray in the same plane and passing over the same streamwise location. Al- 
though this approximation sounds crude at first exposure, it is probably 
satisfactory for most boundary- layer applications, since the contribution to 
q r depends upon the cosine of the angle between the radiant source and the 
normal to the surface and decreases rapidly with distance. Use of the ap- 
proximation affords great simplification in boundary- layer problems since 
q r then depends only on the state of the boundary layer at the local stream- 
wise station. On the basis of these considerations, the one-dimensional 
approximation has been used extensively in boundary- layer studies and seems 
appropriate for use in the present study. 

Absorption has been shown to be important at high reentry velocities 
(above 45,000 feet per second) in an air boundary layer (e.g.. Ref. E-3) . 
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Furthermore, it can be significant at much lower velocities when foreign spe- 
cies are injected by ablation of the surface material. Hence, for generality, 
it is necessary to consider absorption as well as emission. The neglect of 

scattering, on the other hand, is a good assumption under reentry conditions 

E— 2 

as long as solid particles or droplets are not present in the flow field. 

For local thermodynamic equilibrium, the equation of radiative transfer 
E— 1 

can be expressed as 


PH v 61 


I - B (T) 

V V 


(E-l) 


where p is the density, x is the absorption coefficient for the multi- 
component gas, is the specific intensity, and B (T) is Planck's func- 

tion defined by 


B v (T) 


2h 

c 


, 2 (e hvAT _ 1} 


(E-2) 


where h is Planck's constant, c is the velocity of light, v is the fre- 
quency, k is Boltzmann's constant, and T is temperature. The flux of 
radiative energy across the surface da in the frequency interval dv is 
obtained by integrating 1^ cos^duu over all solid angles diu = sin0d$d<p 



(@*cp) cos0sin0dcpd0 


(E-3) 


where 0 is the angle between n (the normal to the surface da) and i 
(the direction of incident radiation) and cp is the angle between a refer- 
ence line in the da plane and the trace of 7 on the da plane (see 
sketch) . 
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In order to allow the incident flux at the edge of the boundary layer to be 
a function of cp and 6 it is expedient to consider the radiant flux in 
the angular interval d Q, q r ^ 0 , obtained by differentiating Eq. (E-3) with 
respect to 0. 


It is also convenient for the purpose of evaluating boundary conditions 
to split the net flux qr v ^ 0 into the contribution q+ v Q in the direction 
of the normal unit vector n and the contribution q“ Q in the opposite 
direction. Then 



(E-4) 


where 


2tt 


v,0 


I V (^# ( P) cos0sin0dcp 0 < 0 < n/2 


2tt 


r v,e 


l v (0,cp) cos0sin0dcp n/2 < 0 < n 


> (E — 5) 


Substituting Eqs. (E-5) into Eq. (E-l) and multiplying by 2nsin0 yields 


d < 4 

- d E3ie = " 2TTSineB^ (T) 0 < e < n/2 

dq r C 

i v.e i = r y.e 


> (E-6) 


pH v dX COS0 


COS0 


2rrsin0B^(T) n/2 < 0 < n J 


Considering the sign convention of the following sketch such that the heat 


boundary- layer edge 


y 


y 



wall 
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flux toward the surface is positive 


dy = - cos0di 


Equations (E-6) can then be written 




P% dy 


d q; 


V,9 


v r e 

COS0 


2nsindB v (T) 


p% dy 


: y,e _ r v, 0 

cos 0 


2nsin0B (T) 
v 7 


0 < 0 < tt/2 


n/2 < 0 < rr J 


Defining an optical thickness t 




P* v dy 


Equations (E-8) become 


dqj <4 

r v,e _ r v,e 

dT COS0 


dq“ q r 

V, 0 _ V, 6 


dT 


cos 0 


- 2rrsin0B^ (T) 


2nsin0B (T) 
v 7 


0 < 0 < n/2 


tt/2 < 0 < TT J 


Solution of Eqs. (E-10) yields 


q + e -t/cos0 

. r v,e 


T 

-l v 


f v 

B (T)sin0e“ t/coS0 dt 0 < 6 < tt/2 




= 2n 


B (T)sin0e t/cos6 dt n / 2 < 0 < tt 


(E-7) 


(E-8) 


(E-9) 


(E-10) 


(E-ll) 


(E-12) 
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where the subscript e refers to the boundary-layer edge. 

It is necessary to evaluate the boundary conditions on the left-hand- 
sides of Eqs. (E-ll) and (E-12) . First, consider the edge boundary condition 
for Eq. (E-12). The first step is to evaluate this equation at the wall. 



v,0 


+ 2tt / B (T) sin0e t//cos0 dt 


n t /cos 6 

e tt/2 < G < n 


(E-13) 


where the subscript w refers to the wall. Assuming diffuse reflection and 
spectral emission from the wall with a transmissivity of zero, the heat flux 
away from the wall is given by E “ 4 



2sin9cos0 




+ rre 


B (T ) 

V V w 

w 


tt/2 < e < tt 

(E-14) 


where a v is the hemispheric surface absorptivity and e v is the hemi- 
spheric surface emissivity and the minus sign arises because of the defini- 
tion of 0. If it is further assumed that the coefficients a Vw and e v for a 
given frequency v depend only upon the nature of the surface and the sur- 
face temperature but not on the radiation field to which the surface is ex- 


(E-15) 


and Eq. (E-14) reduces to 



2sin0cos0 



(1 - e 


v 


) 

w 


+ ne B (T ) 
v v v w' 
w 


tt/2 < 0 < TT 


(E-16) 


The q r which appears in Eq. (E-16) can be expressed in terms of known 
v w 

quantities by evaluating Eq. (E-ll) at the edge. That is 
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-T /cose 


q r = q r eG 

v,0 w v,O e 


+ 2ttJ B v (T)sin0e t / cos0 dt 0 < 0 < rr/2 
0 

(E-17) 


which upon integration yields 


r n//2 -T v /cosd J 1 / 2 ^ v e 

c = q r e 6 d9 + 2n / / - v 

j o V ' 6e j o j o 


B (T)sin0e“ t//cose dtd0 (E-18) 


Exchanging the order of integration in the second term in Eq. (E-18) and 
utilizing the exponential integral, defined by 


E n (t) s / M**' 2e t/U<3 ^ 


(E-19) 


yields the following expression for q + 


w 

- n/ ' 2 -T v / cos0 r v e 

q* e e d0 + 2 tt / B v (T)E 2 (t)dt 




(E-20) 


Introducing Eq. (E-20) into Eq. (E-16) yields the following expression for 
the one-sided heat flux away from the wall evaluated at the wall 


q = 2sin0cos0 

-' e w 


r tt/2 


(1 - e ) 

' V * 

w 


L 0 


-T /cos0 

4 « « 

v,e 


T 

r Ve 1 

1 

+ 2n 1 B v (T)E 2 (t)dt 

+ ne v W 

w 

0 

/ 


n/2 < 0 < n 


(E-21) 
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Substituting Eq. (E-21) into Eq. (E-13) yields the following expression for 
the one-sided heat flux away from the wall evaluated at the boundary- layer 
edge 


/cos 6 


= 2 rre 


v,e 


T 

/ 


B v (T) sin0e" t//cos0 dt 


/cos0 


n/ 2 


- 2sin0cos£(l - e )e 




-T 

e v 


/cos0 1 


30 ' 


/cos0 


4nsin0cos0 (1 - e )e 


B v (T)E 2 (t) dt 


/COS0 

- 2rrsin0cos0e v B^Tje e n/2 < 0 < n (E-22) 

w 

The first through fourth terms in Eq. (E-22) are, respectively, the contri- 
butions from the radiation emitted in the boundary layer at all angles 9 
away from the surface, the radiation incident at the boundary-layer edge 
reflected from the wall, the radiation emitted in the boundary layer at all 
angles 9 toward the surface, and the radiation emitted from the surface. 

... T Ve/cos 6 

The common multiplier e c represents the attenuation* of these types 

of radiation as they pass outward through the boundary layer whereas E 2 (t) 
represents the attenuation of the incident flux as it passes through the 
boundary layer from the edge to the wall. 

The one-sided heat fluxes can now be evaluated. Integration of Eq. 

E-12) yields 


*In that the valid range of 9 is from tt/2 to n, the cos0 is always 
negative . 
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tt/2 


-(t -t )/cos 0 

v v" 

q e d& 

v,0 Q 


TT T 


- 2tt 


B (T) sin0e 


— ( t— T )/cOS0 


dtd0 (E-23) 


tt/2 ^ T 


Substituting for q r ^ q from Eq. (E-22) and inverting the order of integra- 
tion yields 


T TT 

v e C (t -t)/cos0 


q“ = 2 tt / B V (T) / sin0e 

u K/2 


d0dt 


r tt /2 


2(1 - e ) 

v 

w 


L 0 


+ -T /cOS0* Q( 

'Jr fll e v e de 

v,0 4 


T /cOS0 

sin0cos0e v d 0 


J n/2 


4tt (1 - e ) 
w 


r T 


L- 0 


B y (T)E 2 (t) dt 


I sin0cos0e 
tt/2 


T /COS0 


d0 


/ T /COS 0 

sin0cos0e v d0 

n/2 


r v e r -(t-T )/cose 

2n / B v (T) I sinQe d0dt 


(E-24) 


n/2 
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Utilizing the definition of the exponential integrals (Eq. (E-19) ) , Eq. (E-24) 
becomes 


q r = 2 tt / b v (t)e 2 (t v - t) dt 


tt/2 


+ 2(1 - e )E (T ) / q+ e- T v/ COS0 d0 
w J v.e e 


+ 4n(l - e v )E 3 (t v ) / B v (T)E 2 (t)dt 
w / 


+ 2rre B (T ) E 0 ( T ) 
v w v vr 3 ' v 


2rr / B v (T)E 2 (t v - t) dt 


(E-25) 


Integration of Eq. (E-ll) yields for the positive one-sided heat flux 
tt/2 tt/2 t 


► I + T /cose 

C V ■ J \.e u • . , 

0 w 0 0 


■-L I 


B v (T) sinOe 


-(t-T )/cOS0 

v dtdO (E-26) 


Substituting for q + from Eq. (E-17) 


tt/2 


r = q r e 

v j r v,0e 

^ 0 


■(‘ r v -V /cos0 
+ v e v 


t n/2 
d0 + 2 tt I B (T) sin0e 


-(t-T )/cos0 


d0dt 


0 0 


t n/2 
' v r 


- 2tt 


B (T)sin0e 

v 


— (t— T ) /COS0 


d0dt 


(E- 27) 


0 0 
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Again, in terms of exponential integrals 


-(t v -T v )/cose 

q* e e d0 

v '°e 


+ 2n I B v (T)E 2 (t - T v )dt - 2tt I B v (T)E 2 (t - T v ) dt (E-28) 


In accordance with the definition of Eq. (E-4) the net radiant heat 
flux in the frequency interval dv is given by the difference of Eq. (E-28) 
and Eq. (E-25) . Upon combining terms 


2tt B v (T)E 2 (t - T v )dt - 2tt B v (T)E 2 (T v - t) dt 


- 2ne v B v (T w )E 3 (t v ) -toll- e v > E 3<V B v (T)E 2 (t)dt 
w w / 


r / - (T -T )/cOS0 

f + V e V 

+ / q e d0 


- 2d - e v )E (T ) qj • 

w J v,e 

J o e 


-T /cos 0 

v_ 


(E-29) 


The net radiation heat flux, q , is obtained by integrating Eq. (E-29) 
over all frequencies 


q r = / q r dv 

0 v 


(E-30) 
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t 


This integration can be approximated by considering the entire range of v 

to be subdivided into M spectral bands of width 6v . The k and e xi 

m v v w 

are assumed to be frequency independent within a single band but to vary 
from band to band and are thus designated as k m and e m , respectively. 
The choice of the number of bands, M, and the frequency interval of each 
band, can k e ma< ^ e such as to approximate the actual to a given 

degree of accuracy. The q r for this smeared-band model is given by 


M v +6v 


r— t f m m 

= L J q r 

m= 1 v_ v 

m 


dv 


It is convenient to define 


B (T) 
m 


/ 


v +6 

m m 


B (T) dv 
v 


j B v (T) dv 
0 


(E-31) 


(E-32) 


where 


/ B (T) dv = — T 4 

J V tt 

0 


(E-33) 


The net radiant heat flux toward the surface at the nodal point i is then 
given by 



m 


" 4(7(1 " \ )b 3 ( V 


T 4 B m (T)E 2 (t) dt 


tt/ 2 r \ 

l it 


v +6v 
m m 


v.e c 


dv 


( T m ) /cosO 


m 


m 


-T /cose' 


2(1 - e 


■w'W 


)e 


m 


(E-34) 
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■» 


where from Eq. (E-9) T ro is given by 


X =f 


PH m dy 




(E-35) 


At the wall Eq. (E-34) simplifies to 



(E-36) 


since T =0 and E,(t ) =1/2. 
m w 3 m w 

The incident flux at the boundary- layer edge can be a function of fre- 
quency and the angles 0 and cp. Differentiating Eq. (E-3) with respect to 
0 while considering 0 < tt/2 and integrating over the frequency range 6\> m 
yields 


0 < 0 < tt/2 
(E-37) 

The integrations in this term with respect to v and cp can be performed 
a priori for each spectral band. The result is then substituted into Eq. 
(E-34) and (E-35) where the integration with respect to 0 is performed. 

The calculation of q r . proceeds as follows. Given the temperature 
and particle densities across the boundary layer, the T mi and B m (T ± ) ma- 
trices are computed for each spectral band m and for a finite number of 
points across the boundary- layer i (nodal points in the numerical solution 
procedure) . The integrations in Eqs. (E-34) and (E-35) are then performed 
for each spectral band, and the contribution from each spectral band is 
added to yield the q r . 


v +6 v 
m m 


2tt v +6v 
m m 


g* dv 

r v,e e 


= cos0sin0 


m 



I (0,cp)dvdcp 
e 
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APPENDIX F 
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A STUDY OF DIFFERENTIAL VERSUS INTEGRAL PROCEDURES 
EMPLOYING IDENTICAL SPLINE-FIT APPROXIMATIONS 

There are three basic aspects to the integral matrix solution procedure 
described in this report: connected cubics (spline functions) are employed 
to relate the primary dependent variables to the transverse coordinate rj ; 
the differential equations are solved in integral form with a weighting 
function which is unity between neighboring nodal points and zero elsewhere; 
and Newton-Raphson iteration is utilized to solve the resulting set of linear 
and nonlinear equations. The spline-fit approximation of the dependent 
variables was a natural choice since smooth functions were desired to minimize 
the number of nodal points and recent studies have shown spline functions to 
be superior to single higher-order polynomials. Newton-Raphson iteration was 
chosen in order to effect linearized coupling throughout the boundary layer. 
However, the decision of whether or not to integrate and the choice of a 
weighting factor necessitated some study. The results of these studies are 
reported in this appendix. 

F-l 

As pointed out by Dorodnitsyn , solution of the bound ary- layer equations 

in differential form is equivalent to an integral solution using the Dirac 

delta function as the weighting function. The question at hand then resolves 

down to the choice of weighting functions. There are three approaches here: 

the Dirac delta function (differential approach) , a step weighting function 

F — 2 

such as that used by Pallone , or a smooth weighting function across the 

F — 1 

entire boundary layer such as that employed by Dorodnitsyn. “ The primary 
distinction between these three approaches is that the first takes an 
infinitesimally small sample (restricted to the nodal points themselves) , the 
second samples over a portion of the boundary layer, and the third samples 
over the entire boundary layer . 

Parallel developments were made for the first two approaches for non- 
similar incompressible boundary layers, the only distinction being the choice 
of weighting function. The results of this study were inconclusive with 
regard to accuracy, convergence stability, and the number of iterations 
required to achieve convergence. In particular, it was found that the 
accuracy depends more upon the distribution of the nodal points than upon 
the size of the sample. Therefore the choice of weighting function was made 
on the basis of algebraic simplicity which favored step weighting functions 
over Dirac delta functions or smoothly varying functions such as those used 
by Dorodnitsyn. 
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In the remainder of this appendix, Newton-Raphson recurrence formulas 
are presented for the nonsimilar incompressible momentum equation for the 
integral and the differential approaches. Results are then compared for the 
case of an incompressible boundary layer on a flat plate and for a nonsimilar 
boundary layer with an adverse pressure gradient. Finally, the algebraic 
complexity of the two approaches is discussed. 

In the case of an incompressible boundary layer, the energy and species 
conservation equations are not needed, and the transformed momentum equation 
(Eq. (55)) reduces to 


ff" + -f — + P (a!l - f ' s ) 

a H H 
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0 


(F-l) 


The boundary layer at a given streamwise station is divided into N nodal 
points, Tj^ , where i = 1 at the wall and N at the boundary- layer edge, 
in the integral approach, Eq. (F-l) is integrated at constant £ between 
neighboring nodes, ^ and r)^ l 
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(F-2) 


As discussed in the present report, the f^ , f| , f i * and f ' are expanded 
about point i in terms of their ^-derivatives by the use of Taylor series. 
These series are truncated by considering the fV" to be constant between 
T)^ and • Thus, between each i and i+1 the f is represented as 

a quartic, the f' as a cubic , and the f*' as a quadratic, whereas the f ,n is 
considered to vary linearly between each pair of ri-stations, all of these 
functions joining continuously at the nodal points. Utilizing the same pro- 
cedure described in the present report for representing the streamwise deriva- 
tives and evaluating the integrals which appear in Eq. (F-2) , differentiation 
yields the following Newton-Raphson recurrence formulas: 
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where the ERROR is given by the left-hand side of Eq. (F-2) evaluated in the 
previous iteration and the XP^ , XP 2 , XP^ and XP 4 are defined by Eqs. (86) , 
the d Q , d 1 and d 2 by Eqs. (88) or (89) , and the , ZP 2 , ZP 3 and 

ZP 4 by Eq. (94). 

In the differential approach, Eq. (F-l) is solved at each nodal point i 
without integration. Utilizing the same procedure for representing the 
streamwise derivatives, differentiation yields the following Newton-Raphson 
recurrence formulas: 
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Aa R - - ERROR 
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where the ERROR is given by the left-hand side of Eq. (F-l) evaluated in the 
previous iteration. 

Velocity profiles for incompressible flow over a flat plate as reported 
by Howarth F- ^ and as obtained with the above-described integral and matrix 
procedures are presented in Table F-l in terms of the r] defined by Howarth 
(which is times the r\ defined in the present report for this problem). 
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TABLE F-l 


\ 


\ 

VELOCITY PROFILES FOR INCOMPRESSIBLE FLOW OVER A FLAT PLATE 
AFTER HOWARTH COMPARED TO 11 -POINT INTEGRAL-MATRIX 
AND DIFFERENTIAL-MATRIX SOLUTIONS 

Velocity ratio 


r Howarth 

F-3 

Howarth 

Integral-matrix 

Differential- 

0 

0 

0 

0 

.4 

0.13277 

0.1328 

0.1322 

.8 

0.26471 

0.2646 

0.2636 

1.2 

0.39378 

0.3937 

0.3916 

1.6 

0.51676 

0.5166 

0.5139 

2.0 

0.62977 

0.6297 

0.6263 

3.0 

0.84605 

0.8462 

0.8431 

4.0 

0.95552 

0.9553 

0.9564 

5.0 

0.99155 

0.9913 

0.9938 

6.0 

0.99898 

0.9990 

0.9999 

7.0 

0.99992 

0.9999 

1.0000 

0.0 

1.00000 

1.0000 

1.0000 


Identical 11-point nodal distributions were employed in the present calcula- 
tions, namely, j] * 0, 0.2, 0.5, 0.9, 1.4, 2.0, 2.7, 3.5, 4.4, 5.4 and 6.5. 

The values of u/u e reported in Table F-l were then computed from the fl 
and their derivatives by use of the Taylor series expansions (Eqs. (78)). It 
can be seen that the integral-matrix solution agrees with Howarth's results 
to four significant places and that the differential-matrix solution agrees 
within a few tenths of a percent. Corresponding, the wall shear function, 
f^ , was 0.3320 and 0.3308 for the integral and differential procedures, 
respectively, compared to 0.33206 reported by Howarth. 

The problem of linearly retarded flow ( u /u e = 1 - ax with a - 1/8 ) , 

F— 3 ^ F— 4 

first studied by Howarth and later investigated by Smith and Clutter , 

among others, was also considered. Wall shear function, f ^ , as obtained 
with 6-point and 10-point integral and differential solutions are compared in 
Table F-2 to results reported by Smith and clutter. The latter results can 
be considered as precise since they agree closely with those of other investi- 
gators. They obtained this degree of accuracy by the use of small streamwise 
spacing. The present results, on the other hand, were obtained with relatively 
large streamwise spacing. (All stations considered in the present calculations 
are shown whereas only a sampling of the Smith and Clutter results are 
presented. The fact that streamwise spacing is affecting the present results 
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TABLE F-2 


f 


WALL SHEAR FUNCTION ALONG A FLAT PLATE WITH LINEARLY RETARDED FLOW, 

U - 1 - x/8 , AFTER SMITH AND CLUTTER COMPARED TO 10 -POINT AND 6-POINT 
INTEGRAL-MATRIX AND DIFFERENTIAL-MATRIX SOLUTIONS 


Type 

Ref. Differ- 

F-4** ential 

Integral 

Integral 

Differ- 

ential 

Integral 

Integral 

No. Of 
•p-points* 

10 

10 

10 

6 

6 

6 

Normalized 
distance , 

X 

Shear 

function 

evaluated 

at the wall, f M 
* w 


0 

.4696 .4678 

.4695 

.4695 




.01 

.4645 

.4663 

.4663 

.4656 

.4660 

.4660 

.02 



.4635 

.4629 

.4633 

.4633 

.04 

.4555 

.4573 


.4579 

.4582 


.08 

.4451 

.4469 

.4471 

.4472 

.4473 

.4468 

.12 

.4343 

.4361 


.4362 

.4362 


.16 

.4230 

.4249 

.4242 

.4249 

.4248 

.4237 

.20 

.4114 .4114 

.4133 


.4133 

.4130 


.28 

.3866 

. 3886 

.3872 

.3887 

. 3879 

.3862 

.36 

.3601 

.3622 


.3623 

.3610 


.44 

.3315 

.3338 

.3309 

.3340 

.3320 

.3289 

.52 

.3006 

.3030 

.3002 

.3034 

.3002 

.2972 

.60 

.2667 

.2692 

.2667 

.2698 

.2651 

.2623 

.68 

.2288 

.2316 

.2291 

. 2323 

.2251 

.2225 

.7 6 

.1852 

.1882 

.1857 

.1892 

.1778 

.1752 

.80 

.1653 .1605 

.1638 

.1612 

.1649 

.1504 

.1478 

.84 

.1329 

.1365 

.1337 

.1378 

. 1186 

.1160 

.88 

.1100 .1000 

.1042 

.1010 

.1059 

.0809 

.0784 

.92 

.0736 .0506 

.0582 

.0519 

.0615 

.0350 

.0330 

.94 

non*** 


non 

.0162 

non 

.0092 

.96 

0 (extrapj 

non 


non 


-.0132 

.98 






-.0351 

1.00 






-.0549 

1.02 






-.0754 

1.04 






-.0929 

1.06 






-.1135 

* 

10-point ( rj 88 0, 0.2, 0.5, 0 
6-point (rj - 0, 0.8, 1.8, 3. 

.9, 1.4, 
0, 4.4, 6 

2.0, 2.7, 
.0) . 

3.5, 4.4, 

5.4) , 



** 

Smith and Clutter results shown for representative streamwise stations only. 
***nonconvergent 
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can be seen by comparing the two 10-point integral solutions and the two 
6-point integral solutions.) Looking first at the 10-point solutions, the 
integral solution is slightly better except in the vicinity of x = . 20 
where the differential approach agrees with that of Smith and Clutter. In 
the 6-point solutions, the integral and differential methods yield comparable 
results until an x of 0.60 or so, after which the differential approach 
yields better accuracy. 

Since the comparison between the integral and differential approaches 
was inconclusive with regard to accuracy and convergence considerations, it 
was appropriate to consider algebraic complexity. Comparison of Eqs . (F-3) 

and (F-4) indicates that the differential approach is algebraically simpler 
for the incompressible nonsimilar boundary layer. However, in the multi- 
component boundary layer the situation is strikingly reversed. All of the 
complexities introduced by multicomponent thermodynamic and transport proper- 
ties in the energy equation and elemental species equations appear in flux 
divergence terms (see Eqs. (56) and (58)). Integration by a unity weighting 
factor thus eliminates a derivative with respect to q so trial it is not 
necessary to evaluate flux derivatives in the evaluation of the ERRORS. The 
fluxes, in turn, contain derivatives with respect to q (see, e.g., the q* 
given by Eq. (50)). Thus, complete Newton-Raphson iteration requires evalua- 
tion of second derivatives of state functions with respect to h and K-^ 
as well as first derivatives (see Eqs. (129) and (132)). The differential 
approach, on the other hand, would require the evaluation of second and third 
derivatives of state functions with respect to h and for full Newton- 

Raphson iteration. Thus, integration between i and i-1 is by far the most 
desirable approach from the standpoint of algebraic simplicity. 
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